perm filename V2D.IN[1,DEK] blob sn#285515 filedate 1977-05-28 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00015 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	folio 142 galley 1
C00014 00003	folio 144 galley 2
C00025 00004	folio 147 galley 3
C00041 00005	folio 151 galley 4
C00054 00006	folio 156 galley 5
C00069 00007	folio 161 galley 6
C00080 00008	folio 165 galley 7
C00093 00009	folio 173 galley 8
C00106 00010	folio 176 galley 9
C00120 00011	folio 178 galley 10
C00132 00012	folio 182 galley 11
C00145 00013	folio 186 galley 12
C00157 00014	folio 190 galley 13
C00173 00015	folio 194 galley 14
C00185 ENDMK
C⊗;
folio 142 galley 1
    0  {U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addison-W
    1  esley)|9f. 142|9ch. 3|9g. 1c|'{A6}!|9|7In principle, 
    7  any of these other random quantities may be obtained 
   16  from the uniform deviates |εU|β0, U|β1, U|β2,|4.|4.|4.|4. 
   23  |πPeople have devised a number of important ``random 
   31  tricks'' which may be used to perform these manipulations 
   40  e∃ciently on a computer, and a study of these 
   49  techniques also gives some insight into the proper 
   57  use of random numbers in any Monte Carlo application.|'
   66  !|9|7It is conceivable that someday someone will 
   73  invent a random-number generator which produces 
   79  some of these other random quantities |εdirectly, 
   86  |πinstead of getting them indirectly via the 
   93  uniform distribution. Except for the ``random 
   99  bit'' generator described in Section 3.2.2, no 
  106  direct methods have so far proved to be practical.|'
  115  !|9|7The discussion in the following section 
  121  assumes the existence of a random sequence of 
  129  uniformly distributed real numbers between zero 
  135  and one. A new uniform deviate is generated whenever 
  144  we need it. These numbers are normally represented 
  152  in a computer word with the decimal point assumed 
  161  at the left.|'{A18}|∨3|∨.|∨4|∨.|∨1|∨.|9|∨N|∨u|∨n|∨e|∨r|∨i|∨c
  164  |∨a|∨l|4|4|∨D|∨i|∨s|∨t|∨r|∨i|∨b|∨u|∨t|∨i|∨o|∨n|∨s|'
  165  {A6}This section summarizes the best techniques 
  171  known for producing numbers from various important 
  178  distributions. Many of the methods were originally 
  185  suggested by John von Neumann in the early 1950's, 
  194  and these have been gradually improved upon by 
  202  other people, notably George Marsaglia, J. H. 
  209  Ahrens, and U. Dieter.|'{A12}|≡A|≡.|9|≡R|≡a|≡n|≡d|≡o|≡m 
  214  |≡c|≡h|≡o|≡i|≡c|≡e|≡s |≡f|≡r|≡o|≡m |≡a |≡_|≡n|≡i|≡t|≡e 
  218  |≡s|≡e|≡t|≡.|9The simplest and most common type 
  224  of distribution required in practice is a random 
  232  |εinteger. |πAn integer between 0 and 7 can be 
  241  extracted from three bits of |εU |πon a binary 
  250  computer; in such a case, these bits should be 
  259  extracted from the |εmost signi⊂cant (|πleft-hand) 
  265  part of the computer word, since in many random-number 
  274  generators the least signi_cant bits are not 
  281  su∃ciently random. (See the discussion of this 
  288  in Section 3.2.1.1.)|'!|9|7In general, to get 
  295  a random integer |εX |πbetween 0 and |εk|4α_↓|41, 
  303  |πwe can |εmultiply |πby |εk, |πand let |εX|4α=↓|4|"lkU|"L. 
  311  |πOn |¬m|¬i|¬x, we would write|'{A9}|¬l|¬d|¬a!|¬u|;
  317  |¬m|¬u|¬l!|¬k|;{A9}and after these two instructions 
  323  have been executed the desired integer will appear 
  331  in register A. If a random integer between 1 
  340  and |εk |πis desired, we add one to this result. 
  350  [The instruction ``|¬i|¬n|¬c|¬a |¬1'' would follow 
  356  (1).]|'!|9|7This method gives each integer with 
  363  nearly equal probability. (There is a slight 
  370  error because the computer word size is _nite; 
  378  see exercise 2. The error is quite negligible 
  386  if |εk |πis small, for example, if |εk/m|4|π|¬W|41/10000.) 
  394  In a more general situation we might want to 
  403  give di=erent weights to di=erent integers. Suppose 
  410  that the value |εX|4α=↓|4x|β1 |πis to be obtained 
  418  with probability |εp|β1, X|4α=↓|4x|β2 |πwith 
  423  probability |εp|β2,|4.|4.|4.|4, |πand |εX|4α=↓|4x|βk 
  427  |πwith probability |εp|βk. |πWe can generate 
  433  a uniform number |εU |πand let|'{A9}|h|ε|∂X|4α=↓|4|5|6|∂x|βk
  439  ,!!|∂|πif!!|εp|β1|4α+↓|4p|β2|4α+↓|4.|4.|4.|4α+↓|4p|βk|βα_↓|β
  439  1|4|¬Q|4U|4|¬Q|41.|E|n|;|L|Lx|β1,|L|πif!!|ε0|4|¬E|4|εU|4|¬W|
  440  4p|β1;>{A2}{A2}|L{A6}|εX|4α=↓|4{B6}|5|6|Lx|β2,|L|πif!!|εp|β1
  441  |4|¬E|4U|4|¬W|4p|β1|4α+↓|4p|β2;>|L|L|¬O|4|¬O|4|¬O>
  443  |L|L|εx|βk,|L|πif!!|εp|β1|4α+↓|4p|β2|4α+↓|4|¬O|4|¬O|4|¬O|4α+
  443  ↓|4p|βk|βα_↓|β1|4|¬E|4U|4|¬W|41.>{A9}|π(Note 
  445  that |εp|β1|4α+↓|4p|β2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βk|4α=↓|4
  446  1.)|'!|9|7There is a ``best possible'' way to 
  454  do the comparisons of |εU |πagainst various values 
  462  of |εp|β1|4α+↓|4p|β2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βs, 
  464  |πas implied in (2); this situation is discussed 
  472  in Section 2.3.4.5. Special cases can be handled 
  480  by more e∃cient methods; for example, to obtain 
  488  one of the eleven values 2, 3,|4.|4.|4.|4, 12 
  496  with the respective probabilities |f1|d336|), 
  501  |f2|d336|),|4.|4.|4.|4, |f6|d336|),|4.|4.|4.|4, 
  503  |f2|d336|), |f1|d336|), we could compute two 
  509  independent random integers between 1 and 6 and 
  517  add them together.|'!|9|7However, none of the 
  524  above techniques is really the fastest way to 
  532  select |εx|β1,|4.|4.|4.|4,|4x|βk |πwith the correct 
  537  probabilities. A method which is considerably 
  543  more e∃cient in most cases, at a slight increase 
  552  in storage space, is explained below in connection 
  560  with the ``rectangle-wedge-tail'' method and 
  565  further illustrated in exercises 20 and 21. See 
  573  also exercise 28.|'{A12}|≡B|≡.|9|≡G|≡e|≡n|≡e|≡r|≡a|≡l 
  577  |≡m|≡e|≡t|≡h|≡o|≡d|≡s |≡f|≡o|≡r |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o|≡u|
  579  ≡s |≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡i|≡o|≡n|≡s|≡.|9|4The 
  581  most general real-valued distribution may be 
  587  expressed in terms of its ``distribution function'' 
  594  |εF(x), |πwhich speci_es the probability that 
  600  a random quantity |εX |πwill be less than or 
  609  equal to |εx;|'{A6}{A3}F(x)|4α=↓|4|πprobability(|εX|4|¬E|4x)
  612  .|E|;(3)|?{A9}|πThis function always increases 
  618  monotonically from zero to one:|'{A9}|εF(x|β1)|4|¬E|4F(x|β2)
  623  ,!!|πif!!|εx|β1|4|¬E|4x|β2;!!F(|→α_↓|¬X)|4α=↓|40,!!F(|→α+↓|¬
  623  X)|4α=↓|41.|E|;(4)|?{A9}|πExamples of distribution 
  628  functions are given in Section 3.3.1, Fig. 3. 
  636  If |εF(x) |πis continuous and strictly increasing 
  643  (so that |εF(x|β1)|4|¬W|4F(x|β2) |πwhen |εx|β1|4|¬W|4x|β2), 
  648  |πit takes on all values between zero and one, 
  657  and there is an |εinverse function F|gα_↓|g1(y) 
  664  |πsuch that, for 0|4|¬W|4|εy|4|¬W|41,|'{A9}y|4α=↓|4F(x)!!|πi
  668  f|6and|6only|6if!!|εx|4α=↓|4F|gα_↓|g1(y).|E|;
  669  (5)|?{A9}|πA general way to compute a random 
  677  quantity |εX |πwith the continuous, strictly 
  683  increasing distribution |εF(x) |πis to set|'{A9}|εX|4α=↓|4F|
  689  gα_↓|g1(U).|E|;(6)|?{A9}|πFor the probability 
  694  that |εX|4|¬E|4x |πis the probability that |εF|gα_↓|g1(U)|4|
  700  ¬W|4x, |πi.e., the probability that |εU|4|¬W|4F(x), 
  706  |πand this is |εF(x).|'!|9|7|πThe problem now 
  713  reduces to one of numerical analysis, namely 
  720  to _nd good methods for evaluating |εF|gα_↓|g1(U) 
  727  |πto the desired accuracy. Numerical analysis 
  733  lies outsi{U0}{H10L12M29}W58320#Computer Programming|9(Knuth
folio 144 galley 2
  735  /Addison-Wesley)|9f. 144|9ch. 3|9g. 2c|'{A6}!|9|7In 
  740  the _rst place, if |εX|β1 |πis a random variable 
  749  having the distribution |εF|β1(x) |πand if |εX|β2 
  756  |πis an independent random variable with the 
  763  distribution |εF|β2(x), |πthen|'{A9}|∂max(|εX|β1,|4X|β2)!!|∂
  766  |πhas|6the|6distribution!!|εF|β1(x)F|β2(x),|h(x)|4α_↓|4F|β1(
  766  x)F|β2(x).|n|;{A4}|π|Lmin(|εX|β1,|4X|β2)|L|πhas|6the|6distri
  767  bution!!|εF|β1(x)|4α+↓|4F|β2(x)|4α_↓|4F|β1(x)F|β2(x).>
  768  {B20}(7)|E|?{A20}{A9}|π(See exercise 4.) For 
  773  example, the uniform deviate |εU |πhas the distribution 
  781  |εF(x)|4α=↓|4x, |πfor 0|4|¬W|4x|4|¬E|41; |πif 
  785  |εU|β1, U|β2,|4.|4.|4.|4,|4U|βt |πare independent 
  789  uniform deviates, then max(|εU|β1, U|β2,|4.|4.|4.|4,|4U|βt) 
  794  |πhas the distribution function |εF(x)|4α=↓|4x|gt, 
  799  0|4|¬W|4x|4|¬E|41. |πThis is the basis of the 
  806  ``maximum of |εt'' |πtest given in Section 3.3.2. 
  814  Note that the inverse function in this case is 
  823  |εF|gα_↓|g1(y)|4α=↓|4{H11}|¬H{H10}|v4y|). |πIn 
  825  the special case |εtα=↓|π2, we see therefore 
  832  that the two formulas|'{A9}|εX|4α=↓|4{H11}|¬H{H10}U|)!!|πand
  836  !!|εX|4α=↓|4|πmax(|εU|β1,|4U|β2)|E|;(8)|?{A9}|πwill 
  839  give equivalent distributions to the random variable 
  846  |εX, |πalthough this is not obvious at _rst glance. 
  855  We need not take the square root of a uniform 
  865  deviate.|'!|9|7The number of tricks like this 
  872  is endless: any algorithm which employs random 
  879  numbers as input will give a random quantity 
  887  with |εsome |πdistribution as output. The problem 
  894  is to _nd general methods for constructing the 
  902  algorithm, given the distribution function of 
  908  the output. Instead of discussing such methods 
  915  in purely abstract terms, we shall study how 
  923  they can be applied in important cases.|'{A12}|≡C|≡.|9|≡T|≡h
  930  |≡e |≡n|≡o|≡r|≡m|≡a|≡l |≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡i|≡o|≡n|
  932  ≡.|9|4Perhaps the most important nonuniform, 
  937  continuous distribution is the so-called |εnormal 
  943  distribution with mean zero and standard deviation 
  950  one|*/:|'|\{A9}F(x)|4α=↓|4|(1|d2{H11}|¬H{H10}|v22|≤p|)|)|4|↔j
  951  |urx|)α_↓|¬X|)|4e|gα_↓|gt|i2|g/|g2|4dt.|E|;(9)|?
  953  {A9}|πThe signi_cance of this distribution was 
  959  indicated in Section 1.2.10. Note that the inverse 
  967  function |εF|gα_↓|g1 |πis not especially easy 
  973  to compute; but we shall see that several other 
  982  techniques are available.|'{A12}(|*/|↔O|\)|9|εThe 
  986  polar method, |πdue to G. E. P. Box, M. E. Muller, 
  997  and G. Marsaglia.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
 1001  |≡P|9(|εPolar method for normal deviates).|9|πThis 
 1006  algorithm calculates two independent normally 
 1011  distributed variables, |εX|β1 |πand |εX|β2.|'
 1016  {A3}{I1.8H}|π|≡P|≡1|≡.|5|5[Get uniform variables.] 
 1019  Generate two independent random variables, |εU|β1, 
 1025  U|β2, |πuniformly distributed between zero and 
 1031  one. Set |εV|β1|5|¬L|52U|β1|4α_↓|41, V|β2|5|¬L|52U|β2|4α_↓|4
 1034  1. (|πNow |εV|β1 |πand |εV|β2 |πare uniformly 
 1041  distributed between |→α_↓1 and |→α+↓1. On most 
 1048  computers it will be preferable to have |εV|β1 
 1056  |πand |εV|β2 |πrepresented in ⊗oatine-point form 
 1062  at this point.)|'{A3}|≡P|≡2|≡.|9[Compute |εS.] 
 1067  |πSet |εS|5|¬L|5V|=|β1|g2|4α+↓|4V|=|β2|g2.|'{A3}|≡P|≡3|≡.|9[
 1069  |πIs |εS|4|¬R|41?] |πIf |εS|4|¬R|41, |πreturn 
 1074  to step P1. (Steps P1 through P3 are executed 
 1083  1.27 times on the average, with a standard deviation 
 1092  of 0.587; see exercise 6.)|'{A3}|≡P|≡4|≡.|9[Compute 
 1098  |εX|β1, X|β2.] |πSet |εX|β1, X|β2 |πaccording 
 1104  to the following two equations:|'{A9}!!|εX|β1|4α=↓|4V|β1|↔H|
 1109  (|→α_↓2|4|πln|4|εS|d2S|),!!X|β2|4α=↓|4V|β2|↔H|(|→α_↓2|4|πln|
 1109  4|εS|d2S|)|4.|E|;(10)|?{A9}{IC}|πThese are normally 
 1114  distributed variables desired.|'{A12}!|9|7To 
 1118  prove the validity of this method, we use elementary 
 1127  analytic geometry and calculus: If |εS|4|¬W|41 
 1133  |πin step P3, the point in the plane with Cartesian 
 1143  coordinates (|εV|β1,|4V|β2) |πis a |εrandom point 
 1149  uniformly distributed inside the unit circle. 
 1155  |πTransforming to polar coordinates |εV|β1|4α=↓|4R|4|πcos|4|
 1159  ≤], |εV|β2|4α=↓|4R|4|πsin|4|≤], |πwe _nd |εS|4α=↓|4R|g2, 
 1164  X|β1|4α=↓|4{H11}|¬H{H10}|v4|→α_↓2|4|πln|4|εS|)|4|πcos|4|≤], 
 1165  |εX|β2|4α=↓|4{H11}|¬H{H10}|v4|→α_↓2|4|πln|4|εS|)|4|πsin|4|≤]
 1165  . Using also the polar coordinates |εX|β1|4α=↓|4R|¬S|4|πcos|
 1171  4|≤]|¬S, |εX|β2|4α=↓|4R|¬S|4|πsin|4|≤]|¬S, we 
 1174  _nd that |≤]|¬S|4α=↓|4|≤] and |εR|¬S|4α=↓|4{H11}|¬H{H10}|→α_
 1178  ↓2|4|πln|4|εS|). |πIt is clear that |εR|¬S |πand 
 1185  |≤]|¬S are independent, since |εR |πand |≤] are 
 1193  independent inside the unit circle. Also, |≤]|¬S 
 1200  is uniformly distributed between 0 and |ε2|≤p; 
 1207  |πand the probability that |εR|¬S|4|¬E|4r |πis 
 1213  the probability that |→α_↓2|4|πln|4|εS|4|¬E|4r|g2, 
 1217  |πi.e., the probability that |εS|4|¬R|4e|gα_↓|gr|i2|g/|g2. 
 1222  |πThis equals 1|4α_↓|4|εe|gα_↓|gr|i2|g/|g2, |πsince 
 1226  |εS|4α=↓|4R|g2 |πis uniformly distributed between 
 1231  zero and one. The probability that |εR|¬S |πlies 
 1239  between |εr |πand |εr|4α+↓|4dr |πis therefore 
 1245  the derivative of |ε1|4α_↓|4e|gα_↓|gr|i2|g/|g2, 
 1249  |πnamely, |εre|gα_↓|gr|i2|g/|g2|4dr. |πSimilarly, 
 1252  the probability that |≤]|¬S lies between |ε|≤u 
 1259  |πand |ε|≤u|4α+↓|4d|≤u |πis (|ε1/2|≤p)|4d|≤u. 
 1263  |πTherefore the probability that |εX|β1|4|¬E|4x|β1 
 1268  |πand that |εX|β2|4|¬E|4x|β2 |πis|'{A9}{M30}|↔j|ur|)|¬T(|εr,
 1272  |≤u)|3|¬G|3r|4|πcos|4|ε|≤u|¬Ex|β1,r|4|πsin|4|ε|≤u|¬Ex|β2|¬Y|
 1272  )|4|(1|d22|≤p|)|4e|gα_↓|gr|i2|g/|g2|4r|4dr|4d|≤u|4α=↓|4|(1|d
 1272  22|≤p|)|4|↔j|ur|)|¬T(x,y)|3|¬G|3x|¬Ex|β1,y|¬Ex|β2|¬Y|)|4e|gα
 1272  _↓|g(|gx|i2|gα+↓|gy|i2|g)|g/|g2|4dx|4dy|'{A4}{M29}α=↓|4|↔a|↔
 1273  K|(1|d22|≤p|)|4|↔j|urx|β1|)α_↓|¬X|)|4e|gα_↓|gx|i2|g/|g2|4dx|
 1273  ↔s|↔a|↔K|(1|d22|≤p|)|4|↔j|urx|β2|)α_↓|¬X|)|4e|gα_↓|gy|i2|g/|
 1273  g2|4dy|↔s.!(11)|?{A9}|πThis proves |εX|β1 |πand 
 1278  |εX|β2 |πare independent and normally distributed, 
 1284  as desired.|'{A12}|ε(|*/|↔P)|9|\Marsaglia's rectangle-wedge-t
 1287  ail method.|9|4|πIn this method we use the distribution|'
 1295  {A9}|εF(x)|4α=↓|4|↔K|(2|d2|≤p|)|4|↔j|urx|)0|)|4e|gα_↓|gt|i2|
 1295  g/|g2|4dt!!x|4|¬R|40,|E|;(12)|?{A9}|πso that 
 1299  |εF(x) |πgives the distribution of the |εabsolute 
 1306  value |πof a normal deviate. After |εX |πhas 
folio 147 galley 3
 1314  {U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addison-W
 1315  esley)|9f. 147|9ch. 3|9f. 3c|'{A6}!|9|7The rectangle-wedge-t
 1320  ail approach is based on several important general 
 1328  techniques which we shall explore as we develop 
 1336  the algorithm. the _rst key idea is to regard 
 1345  |εF(x) |πas a |εmixture |πof several other functions, 
 1353  namely to write|'{A9}|εF(x)|4α=↓|4p|β1F|β1(x)|4α+↓|4p|β2F|β2
 1356  (x)|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βnF|βn(x),|E|'
 1357  (13)|?{A9}|πwhere |εF|β1, F|β2,|4.|4.|4.|4,|4F|βn 
 1361  |πare appropriate distributions and |εp|β1, p|β2,|4.|4.|4.|4
 1366  ,|4p|βn |πare nonnegative probabilities which 
 1371  sum to 1. If we generate a random variable |εX 
 1381  |πby choosing distribution |εF|βj |πwith probability 
 1387  |εp|βj, |πit is easy to see that |εX |πwill have 
 1397  distribution |εF |πoverall. Some of the distributions 
 1404  |εF|βj(x) |πmay be rather di∃cult to handle, 
 1411  even harder than |εF |πitself, but we can usually 
 1420  arrange things so that the probability |εP|βj 
 1427  |πis very small in this case. Most of the distributions 
 1437  |εF|βj(x) |πwill be quite easy to accommodate, 
 1444  since they will be trivial modi_cations of the 
 1452  uniform distribution. The resulting method yields 
 1458  an extremely e∃cient program, since its |εaverage 
 1465  |πrunning time is very small.|'!|9|7It is easier 
 1473  to understad the method if we work with the |εderivatives 
 1483  |πof the distributions instead of the distributions 
 1490  themselves. Let|'{A6}{A3}|εf|1(x)|4α=↓|4F|¬S(x),!!f|βj(x)|4α
 1492  =↓|4F|=|βj|¬S(x);|;{A9}|πthese are called the 
 1497  |εdensity |πfunctions of the probability distributions. 
 1503  Equation (13) becomes|'{A9}|εf(x)|4α=↓|4p|β1f|β1(x)|4α+↓|4p|
 1506  β2f|β2(x)|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βnf|βn(x).NSS|;
 1507  *?|εf(x)|4α=↓|4p|β1f|β1(x)|4α+↓|4p|β2f|β2(x)|4α+↓|4|¬O|4|¬O|4
 1507  |¬O|4α+↓|4p|βnf|βn(x).|E|;(14)|?{A9}|πEach |εf|βj(x) 
 1511  |πis |→|¬R0, and the total area under the graph 
 1520  of |εf|βj(x) |πso there is a convenient graphical 
 1528  way to display the relation (14): The area under 
 1537  |εf(x) |πis divided into |εn |πparts, with the 
 1545  part corresponding to |εf|βj(x) |πhaving area 
 1551  |εp|βj. |πSee Fig. 9, which illustrates the situation 
 1559  in the case of interest to us here, with |εf(x)|4α=↓|4F|¬S(x
 1568  )|4α=↓|4{H11}|¬H{H10}|v4(2/|≤p)|)|εe|gα_↓|gx|i2|g/|g2; 
 1569  |πthe area under this curve has been divided 
 1577  into |εn|4α=↓|4|π37 parts. There are 12 rather 
 1584  large rectangels, which represent |εp|β1f|β1(x),|4.|4.|4.|4,
 1588   p|β1|β2f|β1|β2(x); |πthere are 12 skinny little 
 1595  rectangles perched on top of these, which represent 
 1603  |εp|β1|β3f|β1|β3(x),|4.|4.|4.|4, p|β2|β4f|β2|β4(x); 
 1605  |πthere are 12 wedge-shaped pieces, which represent 
 1612  |εp|β2|β5f|β2|β5(x),|4.|4.|4.|4, p|β3|β6f|β3|β6(x); 
 1614  |πand the remaining part |εp|β3|β7f|β3|β7(x) 
 1619  |πis the ``tail,'' namely the entire graph of 
 1627  |εf(x) |πfor |εx|4|¬R|43.|'{A6}{H9L11}|π|≡F|≡i|≡g|≡.|7|≡9|≡.
 1630  |9|4The density function divided into 37 parts. 
 1637  The area of each part represents the average 
 1645  number of times a random number with that density 
 1654  is to be computed.|;{A6}{H10L12}!|9|7The rectangular 
 1660  parts |εf|β1(x),|4.|4.|4.|4,|4f|β1|β2(x) |πrepresent 
 1663  |εuniform distributions. |πFor example, |εf|β3(x) 
 1668  |πrepresents a random variable uniformly distributed 
 1674  between |f14Fd32*?|) and |f3|d34|). We want to 
 1681  make it possible to determine the probabilities 
 1688  |εp|β1, p|β2,|4.|4.|4.|4, p|β1|β2 |πe∃ciently 
 1692  as the normal deviate is being generated, so 
 1700  (assuming we have a binary computer#the procedure 
 1707  for a decimal computer is similar) let us take 
 1716  each of these to be multiples of, say, |f1|d3256|). 
 1725  This means we want the area under |εp|βjf|βj(x) 
 1733  |πto be a multiple of |f1|d3256|), so the altitude 
 1742  of |εp|βjf|βj(x) |πis taken to be |"l64|εf(j/4)|"L/64, 
 1749  |πfor |ε1|4|¬E|4j|4|¬E|412. |πThe following table 
 1754  gives the resulting probabilities:|'{A9}|∂!!|∂!!|∂!!|∂!!|∂!!
 1758  |∂!!|∂!!|∂!!|∂!!|∂!!|∂!!|∂!!|∂!!!!!!!!|9|∂!!|5|∂|E|;
 1759  |>|εp|β1|;p|β2|;p|β3|;p|β4|;p|β5|;p|β6|;p|β7|;
 1767  p|β8|;p|β9|;p|β1|β0|;p|β1|β1|;p|β1|β2|;(p|β1|β3|4α+↓|4|¬O|4|
 1772  ¬O|4|¬O|4α+↓p|β3|β7)|;{A7}(15)|?{B7}>{A2}|>|f49|d3256|)|;
 1777  |f45|d3256|)|;|f38|d3256|)|;|f30|d3256|)|;|f23|d3256|)|;
 1781  |f16|d3256|)|;|f11|d3256|)|;|f6|d3256|)|;|f4|d3256|)|;
 1785  |f2|d3256|)|;|f1|d3256|)|;|f0|d3256|)|;|f31|d3256|).|;
 1789  >{A9}|πIt follows that the uniform distributions 
 1796  |εF|β1,|4.|4.|4.|4,|4F|β1|β2 |πmay be used |εp|β1|4α+↓|4|¬O|
 1800  4|¬O|4|¬O|4α+↓ p|β1|β2|4|¬V|487.9 |πpercent of 
 1804  the time; about 88 percent of the area under 
 1813  |εf(x) |πis taken up by the large rectangles.|'
 1821  !|9|7An important technique has been given by 
 1828  Marsaglia for e∃cient selection of the 13 alternatives 
 1836  represented in (15). Given a random binary integer 
 1844  between 0 and 255, whose binary representation 
 1851  is |εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8, |πwe procees 
 1855  as follows:|'{A9}{M29}|πif!|ε!|20|4|¬E|4b|β1b|β2b|β3b|β4|4|¬
 1857  W|410,|'{A2}|πif!|ε|5|540|4|¬E|4|εb|β1b|β2b|β3b|β4b|β5b|β6|4
 1858  |¬W|452,|'{A2}|πif!208|4|¬E|4|εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b
 1859  |β8|4|¬W|4225,|'{A2}|πif!|ε225|4|¬E|4b|β1b|β2b|β3b|β4b|β5b|β
 1860  6b|β7b|β8,|'{B54}{I13.9L}|πlet!|εX|4α=↓|4A[b|β1b|β2b|β3b|β4]
 1861  |4α+↓|4|f1|d34|)U;|'{A2}|πlet!|εX|4α=↓|4B[b|β1b|β2b|β3b|β4b|
 1862  β5b|β6]|4α+↓|4|f1|d34|)U;|'{A2}|πlet!|εX|4α=↓|4C[b|β1b|β2b|β
 1863  3b|β4b|β5b|β6b|β7b|β8]|4α+↓|4|f1|d34|)U;|'{A2}|πgenerate 
 1865  |εX |πusing one of the distributions |εF|β1|β2,|4.|4.|4.|4, 
 1872  F|β3|β7, |πas described below.|'{A9}{IC}{M29}|πHere 
 1877  |εA, B, |πand |εC |πare auxiliary tables which 
 1885  appear in Table 1, and |εU |πis a uniform deviate 
 1895  between zero and one.|'{A9}{H9L11M29}|∂!!!!|6|6|∂!!!!!|∂!!!!
 1899  !|9|∂!!!!!|9|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂|E|;{H8L10}|∨T|∨a|∨b|∨
 1900  l|∨e|4|4|∨1|;{A6}{H9L11}|>|εA[0]|4α=↓|40|'|9A[5]|4α=↓|4|f1|d
 1903  32|)|'|9B[40]|4α=↓|4|f1|d34|)|'|9B[46]|4α=↓|4|f3|d34|)|'
 1906  |9C[208]|4α=↓|40|'|9C[214]|4α=↓|41|'|9C[220]|4α=↓|4|f7|d34|)
 1908  |'>|>A[1]|4α=↓|40|'|9A[6]|4α=↓|4|f1|d32|)|'|9B[41]|4α=↓|4|f1
 1913  |d34|)|'|9B[47]|4α=↓|41|'|9C[209]|4α=↓|4|f1|d34|)|'
 1916  |9C[215]>|>A[2]|4α=↓|40|'|9A[7]|4α=↓|4|f3|d34|)|'
 1920  |9B[42]|4α=↓|4|f1|d34|)|'|9B[48]|4α=↓|4|f3|d32|)|'
 1922  |9C[210]|4α=↓|4|f1|d32|)|'|9C[216]|4α=↓|41|'|9C[222]|4α=↓|4|
 1924  f9|d34|)|'>|>A[3]|4α=↓|4|f1|d34|)|'|9A[8]|4α=↓|41|'
 1929  |9B[43]|4α=↓|4|f1|d32|)|'|9B[49]|4α=↓|4|f3|d32|)|'
 1931  |9C[211]|4α=↓|4|f1|d32|)|'|9C[217]|4α=↓|4|f3|d32|)|'
 1933  |9C[223]|4α=↓|4|f9|d34|)|'>|>A[4]|4α=↓|4|f1|d34|)|'
 1937  |9A[9]|4α=↓|4|f5|d34|)|'|9B[44]|4α=↓|4|f3|d34|)|'
 1939  |9B[50]|4α=↓|4|f7|d34|)|'|9C[212]|4α=↓|4|f3|d34|)|'
 1941  |9C[218]|4α=↓|4|f3|d32|)|'|9C[224]|4α=↓|4|f5|d32|)|'
 1943  >|>|;|;|9B[45]|4α=↓|4|f3|d34|)|'|9B[51]|4α=↓|42|'
 1949  |9C[213]|4α=↓|4|f3|d34|)|'|9C[219]|4α=↓|4|f3|d32|)|'
 1951  >{B2}|J#>{A9}{H10L12M29}|πHere |εA, B, |πand 
 1957  |εC |πare auxiliary tables which appear in Table 
 1965  1, and |εU |πis a uniform deviate between zero 
 1974  and one.|'!|9|7To see why this procedure works, 
 1982  consider for example the case |εj|4α=↓|44. F|β4(x) 
 1989  |πis the uniform distribution between |f3|d34|) 
 1995  and 1, and |εX|4α=↓|4|f3|d34|)|4α+↓|4|f1|d34|)U 
 1999  |πhas this distribution. The probability that 
 2005  this distribution is used in the above method 
 2013  is |f1|d316|) times the number of appearances 
 2020  of ``|f3|d34|)'' in the |εA |πtable, plus |f1|d364|) 
 2028  times the number of its appearances in |εB, |πplus 
 2037  |f1|d3256|) times the number in |εC; |πthis equals 
 2045  |f1|d316|)|4α+↓|4|f3|d364|)|4α+↓|4|f2|d3256|)|4α=↓|4|f30|d32
 2045  56|)|4α=↓|4|εp|β4, |πso |εF|β4(x) |πis selected 
 2050  with the right probability. Of course, it would 
 2058  be even faster to use the following method:|'
 2066  {A9}{M32}``if!!|ε0|4|¬E|4b|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8|4|
 2066  ¬W|4225,!!|πlet!!|εX|4α=↓|4T[b|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β
 2066  8]|4α+↓|4|f1|d34|)U''|'{A9}{M29}|πin place of 
 2070  the _rst three tests in (16). Here |εT |πwould 
 2079  be a table with 225 entries, 49 of which are 
 2089  ``0'', 45 of which are ``|f1|d34|)'', 38 are 
 2097  ``|f1|d32|)'', etc. By comparison, the previous 
 2103  method uses only 39 table entries; and it is 
 2112  only a little slower, since the test ``|εb|β1b|β2b|β3b|β4b|β
 2119  5b|β6|4|¬W|452?'' |πneeds to be made only three-eighths 
 2126  of the time, and the test ``|εb|β1b|β2b|β3b|β4b|β5b|β6b|β6b|
 2132  β8|4|¬W|4225?'' needs to be made only three-sixteenths 
 2139  of the time. The savings of storage space in 
 2148  a method like (16) is even more signi_cant in 
 2157  the more general situation when there are more 
 2165  than 256 possibilities.|'{A6}|π|≡F|≡i|≡g|≡.|7|≡1|≡0|≡.|9|4Fr
 2168  equency functions for which Algorithm L may be 
 2176  used to generate random numbers.|;{A6}{H10L12}!|9|7|πThe 
 2182  reader should pause at this point to study the 
 2191  above method carefully, until he understands 
 2197  why method (16) selects the distributions |εF|βj 
 2204  |πwith the probability |εp|βj, |πfor 1|4|¬E|4j|4|¬E|412.|'
 2210  !|9|7The distributions |εF|β1|β3(x),|4.|4.|4.|4,|4F|β2|β4(x)
 2212   |πare also uniform distributions, and they may 
 2220  be handled similarly. They have quite small probability, 
 2228  however (some number between 0 and |f1|d3256|)-, 
 2235  so it is most e∃cient to test for these probabilities 
 2245  in a part of our normal-deviate-generating alg{U0}{H10L12M29
folio 151 galley 4
 2251  }w58320#Computer Programming|9(Knuth/Addison-Wesley)|9f. 
 2253  151|9ch. 3|9g. 4c|'{A6}!|9|7Now we turn to the 
 2261  question of computing random variables from the 
 2268  wedge-shaped distributions |εF|β2|β5(x),|4.|4.|4.|4, 
 2271  F|β3|β6(x). |πTypical cases are shown in Fig. 
 2278  10; when |εx|4|¬W|41, |πthe curved part is concabe 
 2286  downward, and when |εx|4|¬Q|41 |πit is concave 
 2293  upward, but in each case the curved part is reasonably 
 2303  close to a straight line, and it can be enclosed 
 2313  in two parallel lines as shown.|'{A12}!|9|7In 
 2320  order to handle these wedge-shaped distributions, 
 2326  we will rely on yet another general technique, 
 2334  von Neumann's so-called |εrejection method |πfor 
 2340  obtaining a complicated density from another 
 2346  one which ``encloses'' it. The polar method described 
 2354  above is a simple example of such an approach: 
 2363  Steps P1<P3 obtain a random point inside the 
 2371  unit circle by _rst generating a random point 
 2379  inside the unit circle by _rst generating a random 
 2388  point in a larger square, rejecting it and starting 
 2397  over again if the point was outside the circle.|'
 2406  !|9|7The general rejection method is even more 
 2413  powerful than this. To generate a random variable 
 2421  |εX |πwith density |εf, |πlet |εg |πbe another 
 2429  probability density function such that |εf|1(t)|4|¬E|4cg(t) 
 2435  |πfor all |εt, |πwhere |εc |πis a constant. Now 
 2444  generate |εX |πaccording to density |εg, |πand 
 2451  also generate an independent uniform deviate 
 2457  |εU. |πIf |εU|4|¬R|4f|1(X)/cg(X), |πreject |εX 
 2462  |πand start again with another |εX |πand |εU. 
 2470  |πWhen the condition |εU|4|¬W|4f|1(X)/cg(X) |π_nally 
 2475  occurs, the resulting |εX |πwill have density 
 2482  |εf |πas desired. [Prgoof*?e resulting |εX |πwill 
 2489  have density |εf |πas desired. [Proof: |εX|4|¬E|4x 
 2496  |πwill occur with probability |εp(x)|4α=↓|4|↔j|urx|)0|)|4{H1
 2500  2}({H10}g(t)|4dt|4|¬O|4f|1(t)/cg(t)}h12}){H10}|4α+↓|4gp(x), 
 2501  |πwhere |εq|4α=↓|4|↔j|ur1|)0|)|4{H12}({H10}g(t)|4dt|4|¬O|4(1
 2502  |4α_↓|4f|1(t)/cg(t){H12}){H10}|4α=↓|41|4α_↓|41/c 
 2503  |πis the probability of rejection; hence |εp(x)|4α=↓|4|↔j|ur
 2509  x|)0|)|4f(t)|4dt.]|'!|9|7|πThe rejection technique 
 2513  is most e∃cient when |εc |πis small, since there 
 2522  will be |εc |πiterations on the average before 
 2530  a value is accepted. (See exercise 6.) In some 
 2539  cases |εf|1(x)/cg(x) |πis hard to compute, we 
 2546  may know some bounding functions |εr(x)|4|¬E|4f|1(x)/cg(x)|4
 2551  |¬E|4s(x) |πthat are much simpler, and the exact 
 2559  value of |εf(x)/cg(x) |πneed not be calculated 
 2566  unless |εr(x)|4|¬E|4U|4|¬W|4s(x). |πThe following 
 2570  algorithm solves the ``wedge'' problem by developing 
 2577  the rejection method still further.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i
 2582  |≡t|≡h|≡m |≡L (|εNearly linear densities).|9|4|πThis 
 2587  algorithm may be used to generate a random variable 
 2596  |εX |πfor any distribution whose density |ε(x) 
 2603  |πsatis_es the following conditions (cf. Fig. 
 2609  10):|'{A9}{A8}(18)|E|?{B8}|εf|1(x)|4α=↓|40,!!|πfor!!|εx|4|¬W
 2611  |4s!!|πand|6for!!|εx|4|¬Q|4s|4α+↓|4h;|;{A4}a|4α_↓|4b(x|4α_↓|
 2612  4s)/h|4|¬E|4f|1(x)|4|¬E|4b|4α_↓|4b(x|4α_↓|4s)/h,!!|πfor!!|εs
 2612  |4|¬E|4x|4|¬E|4s|4α+↓|4h.|;{A9}{I1.8H}|π|≡L|≡1|≡.|9[Get 
 2614  |εU|4|¬E|4V.] |πGenerate two independent random 
 2619  variables |εU, V, |πuniformly distributed between 
 2625  zero and one. If |εU|4|¬Q|4V, |πexchenge |εU|6|≠l|6V.|'
 2632  {A3}|π|≡L|≡2|≡.|9[|πEasy case?] If |εV|4|¬E|4a/b, 
 2636  |πgo to L4.|'{A3}|≡L|≡3|≡.|9[Try again?] If |εV|4|¬Q|4U|4α+↓
 2642  |4(1/b)|1f|1(s|4α+↓|4hU), |πgo back to step L1. 
 2648  (If |εa/b |πis close to 1, this step of the algorithm 
 2659  will not be necessary very often.)|'{A3}|≡L|≡4|≡.|9[Compute 
 2666  |εX.] |πSet |εX|5|¬L|5s|4α+↓|4hU.|'{IC}{A9}{H9L11}|π|≡F|≡i|≡
 2669  g|≡.|7|≡1|≡1|≡.|9|4Region of ``acceptance'' in 
 2673  Algorithm L.|;{A9}{H10L12}!|9|7To prove that 
 2678  this algorithm is valid, we observe that when 
 2686  step L4 is reached, the point (|εU,|4V) |πis 
 2694  a random point in the area shaded in Fig. 11, 
 2704  namely, 0|4|¬E|ε|4U|4|¬E|4V|4|¬E|4U|4α+↓|4(1/b)|1f|1(s|4α+↓|
 2705  4hU). |πConditions (18) ensure that|'{A9}|ε|(a|d2b|)|4|¬E|4U
 2710  |4α+↓|4|(1|d2b|)|4f|1(s|4α+↓|4hU)|4|¬E|41.|;{A9}|πNow 
 2712  the probability that |εX|4|¬E|4s|4α+↓|4hx, |πfor 
 2717  0|4|¬E|4x|4|¬E|41, |πis the ratio of area to 
 2724  the left of the vertical line |εU|4α=↓|4x |πin 
 2732  Fig. 11 to the total area, namely,|'{A9}|↔j|ur|εx|)0|)|4|(1|
 2739  d2b|)|4f|1(s|4α+↓|4hu)|4du|↔h|↔j|ur1|)0|)|4|(1|d2b|)|4f|1(s|
 2739  4α+↓|4hu)|4du|4α=↓|4|↔j|ursα+↓hx|)s|)|4f|1(v)|4dv;|;
 2740  {A9}|πtherefore |εX |πhas the correct distribution.|'
 2746  {A9}{H9L11}|≡F|≡i|≡g|≡.|7|≡1|≡2|≡.|9|4The ``rectangle-wedge'
 2747  ' algorithm for generating normal deviates.|;
 2753  {A9}{H10L12}!|9|7With appropriate constants |εa|βj, 
 2757  b|βj, s|βj, |πAlgorithm L will take care of the 
 2766  wedge-shaped densities |εf|βj|βα+↓|β2|β4 |πof 
 2770  Fig. 9, for |ε1|4|¬E|4j|4|¬E|412. |πThe _nal 
 2776  distribution, |εF|β3|β7, |πneeds to be treated 
 2782  only about one time in 370; it is used whenever 
 2792  a result |εX|4|¬R|43 |πis to be given. Exercise 
 2800  11 shows that a standard rejection scheme can 
 2808  be used for this ``tail''; hence we are ready 
 2817  to consider the rectangle-wedge-tail procedure 
 2822  in its entirety:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
 2826  |≡M (|εMarsaglia-MacLaren method for normal deviates).|9|4|π
 2831  This algorithm uses several auxiliary tables, 
 2837  constructed as explained in the text, and examples 
 2845  appear in Tables 1 and 2. We assume that a binary 
 2856  computer is being used; the procedure for a decimal 
 2865  computer is similar.|'{A3}{H10L12}{I2.3H}|5|5|≡M|≡1|≡.|9[Get
 2868   |εU.] |πGenerate a uniform random number |εU|4α=↓|4.b|β0b|β
 2875  1b|β2|4.|4.|4.|4b|βm. (|πHere the |εb'|πs are 
 2880  the bits in the binary representation of |εU. 
 2888  |πFor reasonable accuracy, |εm |πshould be at 
 2895  least 24.) Set |ε|≤c|5|¬L|5b|β0. (|πLater, |ε|≤c 
 2901  |πwill be used to determine the sign of the result.)|'
 2911  {A3}|5|5|≡M|≡2|≡.|9[Large rectangle?] If |εb|β1b|β2b|β3b|β4|
 2914  4|¬W|410, |πwhere ``|εb|β1b|β2b|β3b|β4'' |πdenotes 
 2918  the binary integer |ε8b|β1|4α+↓|44b|β2|4α+↓|42b|β3|4α+↓|4b|β
 2921  4, |πset |εX|5|¬L|5A[b|β1b|β2b|β3b|β4]|4α+↓|4.00b|β5b|β6|4.|
 2923  4.|4.|4b|βm |πand go to M10. Otherwise, if |εb|β1b|β2b|β3b|β
 2930  4b|β5b|β6|4|¬W|452, set |εX|5|¬L|5B[b|β1b|β2b|β3b|β4b|β5b|β6
 2932  ]|4α+↓|4.00b|β7b|β8|4.|4.|4.|4b|βm |πand go to 
 2936  M10. Otherwise, if |εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8|4|¬W|4
 2939  225, |πset |εX|5|¬L|5C[b|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8]|4α+
 2941  ↓|4.00b|β9b|β1|β0|4.|4.|4.|4b|βm |πand go to 
 2945  M10.|'{A9}{IC}{H8L10}|∨T|∨a|∨b|∨l|∨e|4|4|∨2|;
 2947  {A6}{H9L11}EXAMPLE|4|4OF|4|4TABLES|4|4USED|4|4WITH|4|4ALGORI
 2947  THM|4|4Mα/↓|;{M27}{B2}|J#>|∂!!!!!|∂!!!!!|∂!!!!!!!|∂!!!!!!|9|
 2949  ∂!!!!!!|∂!!!!!!|9|∂|E|;|>|εj|;S[j]|;P[j]|;Q[j]|;
 2955  D[j]|;E[j]|;>{B2}|J#>|>|5|51|;0|;0.885|;0.881|;
 2964  0.51|;16!|7|;>|>|5|52|;|f1|d34|)|;0.895|;0.885|;
 2972  0.79|;|5|58!|7|;>|>|5|53|;|f1|d32|)|;0.910|;0.897|;
 2980  0.90|;|5|55.33|;>|>|5|54|;|f3|d34|)|;0.929|;0.914|;
 2988  0.98|;|5|54!|7|;>|>|5|55|;1|;0.945|;0.930|;0.99|;
 2997  |5|53.08|;>|>|5|56|;|f5|d34|)|;0.960|;0.947|;
 3004  0.99|;|5|52.44|;>|>|5|57|;|f3|d32|)|;0.971|;0.960|;
 3012  0.98|;|5|52.00|;>|>|5|58|;|f7|d34|)|;0.982|;0.974|;
 3020  0.96|;|5|51.67|;>|>|5|59|;2|;0.987|;0.982|;0.95|;
 3029  1.43|;>|>10|;|f9|d34|)|;0.991|;0.989|;0.93|;|5|51.23|;
 3038  >|>11|;|f5|d32|)|;0.994|;0.992|;0.94|;|5|51.08|;
 3046  >|>12|;|f11|d34|)|;0.997|;0.996|;0.94|;|5|50.95|;
 3054  >|>13|;3|;1.000|;>{B2}|J#>|Ha*?*?*?{U0}{H10L12M29}{I2.3H}W58320
folio 156 galley 5
 3061  #Computer Programming|9(Knuth/Addison-Wesley)|9f.|4156|9ch. 
 3063  3|9g. 5c|'{A6}|5|5|≡M|≡3|≡.|9[Wedge or tail?] 
 3068  Find the |εsmallest |πvalue of |εj, 1|4|¬E|4j|4|¬E|413, 
 3075  |πfor which |ε.b|β1b|β2|4.|4.|4.|4b|βm|4|¬W|4P[j]. 
 3078  |πIf |εj|4α=↓|413, |πgo to step M8.|'{A3}|5|5|≡M|≡4|≡.|9[Ski
 3084  nny rectangle?] If |ε.b|β1b|β2|4.|4.|4.|4b|βm|4|¬W|4Q[j], 
 3088  |πgenerate a new unidorm random number |εU, |πset 
 3096  |εX|5|¬L|5S[j]|4α+↓|4|f1|d34|)U, |πand go to 
 3100  M10.|'{A3}|5|5|≡M|≡5|≡.|9[Get |εU|4|¬E|4V.] |πGenerate 
 3104  two new uniform deviates, |εU, V; |πif |εU|4|¬Q|4V, 
 3112  |πexchange |εU|5|≠m|5V. (|πWe are now performing 
 3118  Algorithm L.) Set  |εX|5|¬L|5S[j]|4α+↓|4|f1|d34|)U.|'
 3123  {A3}|5|5|π|≡M|≡6|≡.|9[|πEasy case?] If |εV|4|¬E|4D[j], 
 3127  |πgo to M10.|'{A3}|5|5|≡M|≡7|≡.|9[Another try?] 
 3132  If |εV|4|¬Q|4U|4α+↓|4E[j] (e|gα_↓|g(|gX|i2|gα_↓|gS|g[|gj|gα+
 3134  ↓|g1|g]|i2|g)|g/|g2|4α_↓|41), |πgo back to step 
 3139  M5; otherwise go to M10. (This step is executed 
 3148  with low probability.)|'{A3}|5|5|≡M|≡8|≡.|9[Get 
 3152  |εU|g2|4α+↓|4V|g2|4|¬W|41.] |πGenerate two new 
 3156  independent uniform deviates, |εU, V; |πset |εW|5|¬L|5U|g2|4
 3162  α+↓|4V|g2. |πIf |εW|4|¬R|41, |πrepeat this step.|'
 3168  {A3}|5|5|≡M|≡9|≡.|9[Compute |εX|4|¬R|43.] |πSet 
 3171  |εT|4|¬L|4{H11}|¬H{H10}|v4(9|4α_↓|42|4|πln|4|εW)/W|). 
 3172  |πSet |εX|4|¬L|4U|4α⊗↓|4T. |πIf |εX|4|¬R|43, 
 3176  |πgo to M10. Otherwise set |εX|4|¬L|4V|4α⊗↓|4T. 
 3182  |πIf |εX|4|¬R|43, |πgo to M10. Otherwise go back 
 3190  to step M8. (The latter will occur about half 
 3199  the time this step is encountered.)|'{A3}|5|5|≡M|≡8|≡.|9[Get
 3205   supertail deviate.] Generate two new independent 
 3212  uniform deviates, |εU, V, |πand set |εX|4|¬L|4{H11}|¬H{H10}|
 3218  v49|4α_↓|42|4|πln|4|εV|).|'{A3}|5|5|π|≡M|≡9|≡.|9[|πReject 
 3220  ?] If |εUX|4|¬R|43, |πgo back to stop M8. (This 
 3229  will occur only about one-twelfth of the time.)|'
 3237  {A3}|5|5|≡M|≡1|≡0|≡.|9[Attach sign.] If |ε|≤c|4α=↓|41, 
 3241  |πset |εX|4|¬L|4|→α_↓X.|'{IC}{A12}|π!|9|7This 
 3244  algorithm is a very pretty example of mathematical 
 3252  theory intimately interwoven with programming 
 3257  ingenuity#it is a _ne illustration of the art 
 3265  of computer programming.|'!|9|7Tables A, B, and 
 3272  C have already been described; the other tables 
 3280  required by Algorithm M are constructed as follows 
 3288  for 1|4|¬E|4|εj|4|¬E|412, |πin terms of the quantities 
 3295  |εa|βj, b|βj, |πand |εp|βj|βα+↓|β2|β4 |πde_ned 
 3300  in exercise 10:|'{A6}|εP[j]|4α=↓|4p|β1|4α+↓|4p|β2|4α+↓|4|¬O|
 3303  4|¬O|4|¬O|4α+↓|4p|β1|β2|4α+↓|4(p|β1|β3|4α+↓|4p|β2|β5)|4α+↓|4
 3303  |¬O|4|¬O|4|¬O|4α+↓|4(p|βj|βα+↓|β1|β2|4α+↓|4p|βj|βα+↓|β2|β4);
 3303  |'P[13]|4α=↓|41; S[j]|4α=↓|4(j|4α_↓|41)/4; Q[j]|4α=↓|4P[j]|4
 3306  α_↓|4p|βj|βα+↓|β2|β4;|'D[j]|4α=↓|4a|βj/b|βj; 
 3308  E[j]|4α=↓|4{H11}|¬H{H10}2/|≤p|)|4|πexp{H12}({H10}|→α_↓(j/4)|
 3308  g2/2)/b|βjp|βj|βα+↓|β2|β4.|E|'(19)|?|πTable 2 
 3312  shows the values for our example to only a few 
 3322  signi_cant _gures, but in a computer program 
 3329  these quantities would all appear with full-word 
 3336  precision. Algorithm M requires a total of 101 
 3344  auxiliary table entries in all.|'!|9|7This method 
 3351  is extremely fast, since 88 percent of the time 
 3360  only the steps M1, M2, and M10 will need to be 
 3371  executed, and the other steps aren't terribly 
 3378  slow either. We have broken the range from 0 
 3387  to 3 into 12 parts in Fig. 9; if it were broken 
 3399  into more parts, say 48 parts, more table entries 
 3408  would be required, but the computation would 
 3415  stay in steps M1, M2, M10 over 97 percent of 
 3425  the time*3 Complete tables for both binary and 
 3433  decimal computers are given in the article by 
 3441  Marsaglia, MacLaren, and Bray, |εCACM|67 (1964), 
 3447  4<10; |πan additional trick, overlapping parts 
 3453  of Tables A, B, C and S, is employed there in 
 3464  order to save storage space. Further re_nements 
 3471  have been developed by Marsaglia, Anathanarayanax, 
 3477  and Paul, |εInf. Proc. Letters |≡4 (1976), |πto 
 3485  appear.|'{A12}|ε|*/(|↔L) |\Forsythe's method.|9|4|πAn 
 3489  amazingly simple technique for generating random 
 3495  deviates with a density of the general exponential 
 3503  form|'{A9}|εf|1(x)|4α=↓|4Ce|gα_↓|gh|g(|gx|g),!!|πfor!!|εa|4|
 3504  ¬E|4x|4|¬W|4b;|6f|1(x)|4α=↓|40|6|πotherwise;|E|;
 3505  (20)|?{A4}|ε0|4|¬E|4h(x)|4|¬E|41!!|πfor!!|εa|4|¬E|4x|4|¬W|4b
 3506  ;|E|;(21)|?{A9}|πwas discovered by John von Neumann 
 3514  and G. E. Forsythe about 1950. The idea is based 
 3524  on the rejection method described earlier, letting 
 3531  |εg(x) |πbe the uniform distribution on [|εa,b): 
 3538  |πWe set |εX|4|¬L|4a|4α+↓|4(b|4α_↓|4a)U, |πwhere 
 3542  |εU |πis a uniform deviate, and then we want 
 3551  to accept |εX |πwith probability |εe|gα_↓|gh|g(|gX|g). 
 3557  |πThe latter operation could be done by testing 
 3565  |εe|gα_↓|gh|g(|gX|g) |πvs. |εV, |πor |εh(X) |πvs. 
 3571  |ε|→α_↓|πln |εV, |πwhen |εV |πis another uniform 
 3578  deviate, but the job can be done without applying 
 3587  any transcendental functions in the following 
 3593  interesting way. Set |εV|β0|4|¬L|4h(X), |πthen 
 3598  generate uniform deviates |εV|β1, V|β2,|4.|4.|4. 
 3603  |πuntil _nding some |εK|4|¬R|41 |πwith |εV|βK|βα_↓|β1|4|¬W|4
 3608  V|βK. |πFor _xed |εV |πand |εk, |πthe probability 
 3616  that |εH(X)|4|¬R|4V|β1|4|¬R|4|¬O|4|¬O|4|¬O|4|¬R|4V|βk 
 3618  |πis |ε1/k*3 |πtimes the probability that max|4(|εV|β1,|4.|4.
 3624  |4.|4,|4V|βk){H10}|4|¬E|4|εh(X), |πnamely |εh(X)|gk/k*3; 
 3627  |πhence the probability that |εK|4α=↓|4k |πis 
 3633  |εh(X)|gk|gα_↓|g1/(k|4α_↓|41)*3|4α_↓|4h(X)|gk/k*3, 
 3634  |πand the probability that |εK |πis odd is|'{A9}|ε|↔k|uc|)k|
 3642  4|πodd|d|εk|¬R1|)|4|↔a|(h(X)|gk|gα_↓|g1|d2(k|4α_↓|41)*3|)|4α_
 3642  ↓|4|(h(X)|gk|d2k*3|)|↔s|4α=↓|4e|gα_↓|gh|g(|gX|g).|;
 3643  {A9}|πTherefore we reject |εX |πand try again 
 3650  if |εK |πis even; we accept |εX |πas a random 
 3660  variable with density (20) if |εK |πis odd. Note 
 3669  that we usually won't have to generate many |εV'|πs 
 3678  in order to determine |εK, |πsince the average 
 3686  value of |εK (|πgiven |εX) |πis |¬K|ur|)|εk|¬R0|)|4|πPr(|εK|
 3692  4|¬Q|4k)|4α=↓|4|¬K|ur|)k|¬R0|)|4h(X)|gk/k*3|4α=↓|4e|gh|g(|gX|
 3692  g)|4|¬E|4e.|'!|9|7|πForsythe realized some years 
 3697  later that this approach leads to an e∃cient 
 3705  method for calculating normal deviates, without 
 3711  the need for any |πauxiliary routines to calculate 
 3719  square roots or logarithms as in Algorithms P 
 3727  and M. His procedure, with an improved choice 
 3735  of intervals [|εa,b) |πdue to J. H. Ahrens and 
 3744  U. Dieter, can be summarized as follows.|'{A12}|≡A|≡l|≡g|≡o|
 3751  ≡r|≡i|≡t|≡h|≡m |≡F (|εForsythe method for normal 
 3757  deviates).|9|4|πThis algorithm generates normal 
 3761  deviates on a binary computer, assuming approximately 
 3768  |εm|4α+↓|41 |πbits of accuracy. |πExtremely large 
 3774  values of |εX |πwhose probability is less than 
 3782  2|ε|gα_↓|gm|gα_↓|g1 |πwill never occur. The algorithm 
 3788  requires a table of values |εd|βj|4α=↓|4a|βj|4α_↓|4a|βj|4α_↓
 3793  |41, |πfor |ε1|4|¬E|4j|4|¬E|4m|4α+↓|41, |πwhere 
 3797  |εa|βj |πis de_ned by the relation|'{A9}|ε|↔H|(2|d2|≤p|)|4|↔
 3803  j|ur|¬X|)a|βj|)|4e|gα_↓|gt|i2|g/|g2|4dt|4α=↓|4|(1|d22|gj|)|4
 3803  .|E|;(23)|?{A9}{I1.8H}|π|≡F|≡1|≡.|9[Get |εU.] 
 3807  |πGenerate a uniform random number |εU|4α=↓|4.b|β0b|β1|4.|4.
 3812  |4.|4b|βm, |πwhere |εb|β0,|4.|4.|4.|4,|4b|βm 
 3815  |πdenote the bits in binary notation. Set |ε|≤c|4|¬L|4b|β0, 
 3823  j|4|¬L|41, |πand |εa|4|¬L|40.|'{A3}|π|≡F|≡2|≡.|9[Find 
 3827  _rst zero |εb|βj.]|9|πIf |εb|βj|4α=↓|41, |πset 
 3832  |εa|4|¬L|4a|4α+↓|4d|βj, j|4|¬L|4j|4α+↓|41, |πand 
 3835  |εrepeat |πthis step. (If |εj|4α=↓|4m|4α+↓|41, 
 3840  |πtreat |εb|βj |πas zero.)|'{A3}|≡F|≡3|≡.|9[Generate 
 3845  candidate.]|9(Now |εa|4α=↓|4a|βj|βα_↓|β1, |πand 
 3848  the current value of |εj |πoccurs with probability 
 3856  |→|¬V1/2|ε|gj. |πWe will generate |εX |πin the 
 3863  range [|εa|βj|βα_↓|β1,|4a|βj), |πwing the rejection 
 3868  method described above, with |εh(x)|4α=↓|4x|g2/2|4α=↓|4y|g2/
 3872  2|4α+↓|4ay |πwhere |εy|4α=↓|4x|4α_↓|4a. |πExercise 
 3876  12 proves that |εh(x)|4|¬E|41 |πas required in 
 3883  (21).) Set |εY|4|¬L|4d|βj |πtimes |ε.b|βj|βα+↓|β1|4.|4.|4.|4
 3887  b|βm |πand |εV|4|¬L|4(|f1|d32|)Y|4α+↓|4a)Y. (|πSince 
 3891  the average value of |εj |πis 2, there will usually 
 3901  be enough signi_cant bits in |ε.b|βj|βα+↓|β1|4.|4.|4.|4b|βm 
 3907  |πto provide decent accuracy. The calculations 
 3913  are readily done in _xed point arithmetic.)|'
 3920  {A3}|≡F|≡4|≡. [Reject*3]|9|4Generate a uniform 
 3924  deviate |εU. |πIf |εV|4|¬W|4U, |πgo on to step 
 3932  F5. Otherwise set |εV |πto a new uniform deviate; 
 3941  and if now |εU|4|¬W|4V (|πi.e., if |εK|πis even 
 3949  in the discussion above), go back to F3, otherwise 
 3958  repeat s{U0}{H10L12M29}|πW58320#Computer Programming|9(Knuth
folio 161 galley 6
 3960  /Addison-Wesley)|9f. 161|9ch. 3|9g. 6c|'{A6}!|9|7Values 
 3965  of |εd|βj |πfor |ε1|4|¬E|4j|4|¬E|447 |πappear 
 3970  in a paper by Ahrens and Dieter, |εMath. comp. 
 3979  |π|≡2|≡7 (1973), 927<937. This paper discusses 
 3985  re_nements of the algorithm which improve its 
 3992  speed at the expense of more tables. Algorithm 
 4000  F is attractive since it is almost as fast as 
 4010  Algorithm M and it is much easier to implement. 
 4019  The average number of uniform deviates required 
 4026  per normal deviate is 2.53947; Brent [(|εACM 
 4033  |π|≡1|≡7 (1974), 704<705] has shown how to reduce 
 4041  this number to 1.37446 at the expense of two 
 4050  subtractions and one division per uniform deviate 
 4057  saved.|'{A12}|ε|*/(|↔C)|9|\Variations of the normal 
 4062  distribution.|9|4|πSo far we have considered 
 4067  the normal distribution with mean zero and standard 
 4075  deviation one. If |εX |πhas this distribution, 
 4082  then|'{A9}|εY|4α=↓|4|≤m|4α+↓|4|≤sX|E|;(24)|?{A9}|πhas 
 4086  the normal distribution with mean |ε|≤m |πand 
 4093  standard deviation |ε|≤s. |πFurthermore, if |εX|β1 
 4099  |πand |εX|β2 |πare independent normal deviates 
 4105  with mean zero and standard deviation one, and 
 4113  if|'{A9}|εY|β1|4α=↓|4|≤m|β1|4α+↓|4|≤s|β1X|β1,!!Y|β2|4α=↓|4|≤
 4114  m|β2|4α+↓|4|≤s|β2(|≤rX|β1|4α+↓|4{H11}|¬H{H10}|v41|4α_↓|4|≤r|
 4114  g2|)X|β2),|E|;(25)|?{A9}|πthen |εY|β1 |πand |εY|β2 
 4120  |πare |εdependent |πrandom variables, normally 
 4125  distributed with menas |ε|≤m|β1, |≤m|β2 |πstandard 
 4131  deviations |ε|≤s|β1, |≤s|β2, |πand with correlation 
 4137  coe∃cient |ε|≤r. (|πFor a generalization to |εn 
 4144  |πvariables, see exercise 13.)|'{A12}|≡D|≡.|9|≡T|≡h|≡e 
 4149  |≡e|≡x|≡p|≡o|≡n|≡e|≡n|≡t|≡i|≡a|≡l |≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡
 4150  t|≡i|≡o|≡n|≡.|9|4After uniform deviates and normal 
 4155  deviates, the next most important random quantity 
 4162  is an |εexponential deviate. |πSuch numbers occur 
 4169  in ``arrival time'' situations; for example, 
 4175  if a radioactive substance emits alpha particles 
 4182  at a rate so that one particle is emitted every 
 4192  |ε|≤m |πseconds on the average, then the time 
 4200  between two successive emissions has the exponential 
 4207  distribution with mean |ε|≤m. |π|πThis distribution 
 4213  is de_ned by the formula|'{A9}|εF(x)|4α=↓|41|4α_↓|4e|gα_↓|gx
 4218  |g/|g|≤m,!!x|4|¬R|40.|E|;(26)|?{A9}{A3}|*/(|↔O)|9|\Logarithm 
 4221  method.|9|4|πClearly, if |εy|4α=↓|4F(x)|4α=↓|41|4α_↓|4e|gα_↓
 4223  |gx|g/|g|≤m, x|4α=↓|4F|gα_↓|g1(y)|4α=↓|4|→α_↓|≤m|πln(|ε1|4α_
 4224  ↓|4y). |πTherefore by Eq. (6), |→α_↓|ε|≤m|πln(|ε1|4α_↓|4U) 
 4230  |πhas the exponential distribution. Since |ε1|4α_↓|4U 
 4236  |πis uniformly distributed when |εU |πis, we 
 4243  conclude that|'{A9}|εX|4α=↓|4|→α_↓|≤m|πln|4|εU|E|;
 4246  (27)|?{A9}|πis exponentially distributed with 
 4251  mean |ε|≤m. (|πThe case |εU|4α=↓|40 |πmust be 
 4258  avoided.)|'{A12}|ε|*/(|↔P)|9|\Random minimization 
 4261  method.|9|4|πWe saw in Algorithm F that there 
 4268  are simple and fast alternatives to calculating 
 4275  the logarithm of a uniform deviate. The following 
 4283  especially e∃cient approach has been developed 
 4289  by G. Marsaglia, M. Sibuya, and J. H. Ahrens.|'
 4298  {A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡S (|εExponential 
 4301  distribution with mean |≤m).|9|4|πThis algorithm 
 4306  produces exponential deviates on a binary computer, 
 4313  using uniform deviates of |εm-|πbit accuracy. 
 4319  The constants|'{A9}|εQ[k]|4α=↓|4|(|πln|42|d21*3|)|4α+↓|4|((ln
 4321  |42)|g2|d22*3|)|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4|((ln|42)|g|εk|d2k*3
 4321  |)|4,!!k|4|¬R|41,|E|;(28)|?{A9}|πshould be precomputed, 
 4326  extending until |εQ[k]|4|¬Q|41|4α_↓|42|g1|gα_↓|gm.|'
 4329  {A3}{I1.7H}|π|≡S|≡1|≡.|9[Get |εU |πand shift.] 
 4333  Generate a uniform random binary fraction |εU|4α=↓|4.b|β1b|β
 4339  2|4.|4.|4.|4b|βm; |πLocate the _rst zero bit 
 4345  |εb|βj, |πand shift o= the leading |εj |πbits, 
 4353  setting |εU|5|¬L|5.b|βj|βα+↓|β1|4.|4.|4.|4b|βm. 
 4355  (|πAs in Algorithm F, the average value of |εj 
 4364  |πis 2.)|'{A3}|≡S|≡2|≡.|9[Immediate acceptance?]|9|4If 
 4368  |εU|4|¬W|4|πln|42, set |εX|5|¬L|5|ε|≤m(j|4|πln|42|4α+↓|4|εU)
 4370   |πand terminate the algorithm. (Note that |εQ[1]|4α=↓|4|πln
 4377  |42.)|'{A3}|≡S|≡3|≡.|9[Minimize.]|9|4Find the 
 4380  least |εk|4|¬R|42 |πsuch that |εU|4|¬W|4Q[k]. 
 4385  |πGenerate |εk |πnew uniform deviates |εU|β1,|4.|4.|4.|4,|4U
 4390  |βk |πand set |εV|5|¬L|5|πmin(|εU|β1,|4.|4.|4.|4,|4U|βk).|'
 4394  {A3}|≡S|≡4|≡.|9[Deliver the answer.]|9|4|πSet 
 4397  |εX|5|¬L|5|≤v(j|4α+↓|4V)|πln|42.|'{A12}{IC}|≡E|≡.|9|≡O|≡t|≡h
 4398  |≡e|≡r |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o|≡u|≡s |≡d|≡i|≡s|≡t|≡r|≡i|≡b|
 4400  ≡u|≡t|≡i|≡o|≡n|≡s|≡.|9|4Let us now consider brie⊗y 
 4405  how to handle some other distributions which 
 4412  arise reasonably often in practice.|'{A12}|ε|\(|↔O)|9|\The 
 4418  gamma distribution |πof order |εa|4|¬Q|40 |πis 
 4424  de_ned by|'{A9}|εF(x)|4α=↓|4|(1|d2|≤)(a)|)|4|↔j|urx|)0|)|4t|
 4426  ga|gα_↓|g1e|gα_↓|gtdt,!!x|4|¬R|40.|E|;(29)|?{A9}|πWhen 
 4429  |εa|4α=↓|41, |πthis is the exponential distribution 
 4435  with mean 1; when |εa|4α=↓|4|f1|d32|), |πit is 
 4442  the distribution of |f1|d32|)|εZ|g2, |πwhere 
 4447  |εZ |πhas the normal distribution (mena 0, variance 
 4455  1). If |εX |πand |εY |πare independent gamma-distributed 
 4463  random variables, of order |εa |πand |εb |πrespectively, 
 4471  then |εX|4α+↓|4Y |πhas the gamma distribution 
 4477  of order |εa|4α+↓|4b. |πThus, for example, the 
 4484  sum of |εk |πindependent exponential deviates 
 4490  with mean 1 has the gamma distribution of order 
 4499  |εk. |πIf the logarithm method (27) is being 
 4507  used to generate these exponential deviates, 
 4513  we need compute only one logarithm: |εX|4|¬L|4|→α_↓|πln(|εU|
 4519  β1|4.|4.|4.|4U|βk), |πwhere |εU|β1,|4.|4.|4.|4,|4U|βk 
 4522  |πare nonzero uniform deviates. This technique 
 4528  handles all integer orders |εa; |πto complete 
 4535  the picture, a suitable method for 0|4|¬W|4a|4|¬W|41 
 4542  |πappears in exercise 16.|'!|9|7The simple logarithm 
 4549  method is much too slow when |εa |πis large, 
 4558  since it requires |"l|εa|"L |πuniform deviates. 
 4564  For large |εa, |πthe following remarkable algorithm 
 4571  due to J. H. Ahrens is reasonably e∃cient, and 
 4580  easy to write in terms of standard subroutines.|'
folio 165 galley 7
 4588  |{U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addison-
 4589  Wesley)|9f. 165|9ch. 3|9g. 7c|'{A6}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|
 4593  ≡m |≡A (|εGamma distribution of order a|4|¬Q|41).|'
 4600  {A3}{I1.8H}|π|≡A|≡1|≡.|9[Generate candidate.] 
 4602  Set |εY|4|¬L|4|πtan(|ε|≤pU), |πwhere |εU |πis 
 4607  a uniform deviate, and set |εX|5|¬L|5{H11}|¬H{H10}|v42a|4α_↓
 4612  |41|)|2Y|4α+↓|4a|4α_↓|41. (|πIn place of tan(|ε|≤pU) 
 4617  |πwe could use a polar method, e.g., |εV|β2/V|β1 
 4625  |πin step P4 of Algorithm P.)|'{A3}|≡A|≡2|≡.|9[Accept?] 
 4632  If |εX|4|¬E|40, |πreturn to A1. Otherwise generate 
 4639  a new uniform deviate |εU, |πand return to A1 
 4648  if |εU|4|¬Q|4(1|4α+↓|4Y|g2)|πexp{H12}({H10}(|εa|4α_↓|41)|πln
 4649  {H12}({H10}X/(a|4α_↓|41){H12}){H10}|4α_↓|4{H11}|¬H{H10}|v42a
 4649  |4α_↓|41|)|2Y{H12}){H10}.|'{A12}{IC}|πThe average 
 4652  number of times step A1 is performed is |→|¬W1.902 
 4661  when |εa|4|¬R|43. |πFor further discussion, proof, 
 4667  and a slightly more complex method which is two 
 4676  to three times faster, see J. H. Ahrens and U. 
 4686  Dieter, |εComputing |π|≡1|≡2 (1974), 223<246. 
 4691  A still faster method, using normal and exponential 
 4699  deviates, has been developed by G. Marsaglia 
 4706  [to appear].|'{A12}|ε|*/(|"P)|9|\The beta distribution 
 4711  |πwith positive parameters |εa |πand |εb |πis 
 4718  de_ned by|'{A9}|εF(x)|4α=↓|4|(|≤)(a|4α+↓|4b)|d2|≤)(a)|≤)(b)|
 4720  )|4|↔j|urx|)0|)|4t|ga|gα_↓|g1(1|4α_↓|4t)|gb|gα_↓|g1|4dt,!!0|
 4720  4|¬E|4x|4|¬E|41.|E|;(30)|?{A9}|πLet |εX|β1 |πand 
 4725  |εX|β2 |πbe independent gamma deviates of order 
 4732  |εa |πand |εb, |πrespectively, and set |εX|5|¬L|5X|β1/(X|β1|
 4738  4α+↓|4X|β2). |πAnother method, useful for small 
 4744  |εa |πand |εb, |πis to set |εY|β1|5|¬L|5U|ur1/a|)1|) 
 4751  |πand |εY|β2|5|¬L|5U|ur1/b|)2|) |πrepeatedly 
 4754  until |εY|β1|4α+↓|4Y|β2|4|¬E|41; |πthen |εX|5|¬L|5Y|β1/(Y|β1
 4757  |4α+↓|4Y|β2). [|πSee M. D. J|=4ohnk, |εMetrika 
 4763  |π|≡8 (1964), 5<15.] |πStill another approach, 
 4769  if |εa |πand |εb |πare integers (not too large), 
 4778  is to set |εX |πto the |εb|πth largest of |εa|4α+↓|4b|4α_↓|4
 4787  1 |πindependent uniform deviates (cf. exercise 
 4793  5<7).|'{A12}|ε|*/(|↔L)|9|\The chi-square distribution 
 4797  |πwith |ε|≤n |πdegrees of freedom (Eq. 3.3.1<22) 
 4804  is obtained by setting |εX|5|¬L|52Y, |πwhere 
 4810  |εY |πhas the gamma distribution of order |ε|≤n/2.|'
 4818  {A12}|ε|*/(|↔C)|9|\The F-distribution (|πvariance-ratio 
 4821  distribution) with |εv|β1 |πand |εv|β2 |πdegrees 
 4827  of freedom, de_ned by|'{A9}|εF(x)|4α=↓|4|(|≤n|ur|≤n|β1/2|)1|
 4831  )|≤n|ur|≤n|β2/2|)2|)|≤){H12}({H10}(|≤n|β1|4α+↓|4|≤n|β2)/2{H1
 4831  2}){H10}|d2|≤)(|≤n|β1/2)|≤)(|≤n|β2/2)|)|4|↔j|urx|)0|)|4t|g|≤
 4831  n|r1|g/|g2|gα_↓|g1(|≤n|β2|4α+↓|4|≤n|β1t)|gα_↓|g|≤n|r1|g/|g2|
 4831  gα_↓|g|≤n|r2|g/|g2|4dt,|E|;(31)|?{A9}|πwhere 
 4834  |εx|4|¬R|40. |πLet |εY|β1 |πand |εY|β2 |πbe independent, 
 4841  having the chi-square distribution with |ε|≤n|β1, 
 4847  |≤n|β2 |πdegrees of freedom, respectively; set 
 4853  |εX|5|¬L|5Y|β1|≤n|β2/Y|β2|≤n|β1. |πOr set |εX|5|¬L|5|ε|≤n|β2
 4856  Y/|≤n|β1(1|4α_↓|4Y), |πwhere |εY |πhas the beta 
 4862  distribution with parameters |ε|≤n|β1/2, |≤n|β2/2.|'
 4867  {A12}|ε|*/(|↔C)|9|\The t-distribution |πwith |ε|≤n 
 4871  |πdegrees of freedom, de_ned by|'{A9}|εF(x)|4α=↓|4|(|≤){H12}
 4876  ({H10}(|≤n|4α+↓|41)/2{H12}){H10}|d2{H11}|¬H{H10}|v2|≤p|≤n|)|
 4876  4|≤)(|≤n/2)|)|4|↔j|urx|)α_↓|¬X|)|4(1|4α+↓|4t|g2/|≤n)|gα_↓|g(
 4876  |g|≤n|gα+↓|g1|g)|g/|g2|4dt.|E|;(32)|?{A9}|πLet 
 4879  |εY|β1 |πbe a nromal deviate (mean 0, variance 
 4887  1) and let |εY|β2 |πbe independent of |εY|β1, 
 4895  |πhaving the chi-square distribution with |ε|≤n 
 4901  |πdegrees of freedom; set |εX|5|¬L|5Y|β1/{H11}|¬H{H10}|v4Y|β
 4905  2/|≤n|).|'{A12}|*/(|↔o)|9|\Random point on n-dimensional 
 4910  sphere with radius one.|9|4|πLet X|β1, X|β2,|4.|4.|4.|4, 
 4916  X|βn |πbe independent normal deviates (mean 0, 
 4923  variance 1); the desired point on the unit sphere 
 4932  is|'{A9}|ε(X|β1/r,|6X|β2/r,|4.|4.|4.|4,|6X|βn/r),!!|πwhere!!
 4933  |εr|4α=↓|4{H11}|¬H{H10}|v4X|=|β1|g2|4α+↓|4X|=|β2|g2|4α+↓|4|¬
 4933  O|4|¬O|4|¬O|4α+↓|4X|=|βn|g2|).|E|;(33)|?{A9}|πNote 
 4936  that if the |εX'|πs are calculated using the 
 4944  polar method, Algorithm P, we compute two independent 
 4952  |εX'|πs each time, and |εX|=|β1|g2|4α+↓|4X|=|β2|g2|4α=↓|4|→α
 4956  _↓2|4|πln|4|εS (|πin the notation of that algorithm); 
 4963  this saves a little of the time needed to evaluate 
 4973  |εr. |πThe validity of this method comes from 
 4981  the fact that the distribution function for the 
 4989  point (|εX|β1,|4.|4.|4.|4, X|βn) |πhas a density 
 4995  which depends only on the distance from the origin, 
 5004  so when it is projected onto the unit sphere 
 5013  it has the uniform distribution. this method 
 5020  was _rst suggested by G. W. Brown, in |εModern 
 5029  Mathematics for the Engineer, |πFirst series, 
 5035  ed. by E. F. Beckenbach (New York: McGraw-Hill, 
 5043  1956), p. 302. To get a random point |εinside 
 5052  |πthe |εn-|πsphere, R. P. Brent suggests taking 
 5059  a point on the surface and multiplying it by 
 5068  |εU|g1|g/|gn.|'|π!|9|7In three dimensions a signi_cantly 
 5074  simpler method can be used, since each individual 
 5082  coordinate is uniformly distributed between |→α_↓1 
 5088  and 1: Find |εV|β1, V|β2, S |πby steps P1<P3 
 5097  of Algorithm P; then the desired point is (|ε|≤aV|β1, 
 5106  |≤aV|β2, 2S-1), |πwhere |ε|≤a|4α=↓|42{H11}|¬H{H10}|v41|4α_↓|
 5109  4S|). [|πRobert E. Knop, |εCACM |π|≡1|≡3 (1970), 
 5116  326.]|'{A12}|≡F|≡.|9|≡I|≡m|≡p|≡o|≡r|≡t|≡a|≡n|≡t 
 5118  |≡i|≡n|≡t|≡e|≡g|≡e|≡r|≡-|≡v|≡a|≡l|≡u|≡e|≡d |≡d|≡i|≡s|≡t|≡r|≡
 5119  i|≡b|≡u|≡t|≡i|≡o|≡n|≡s|≡.|9|4A probability distribution 
 5122  which consists only of integer values can essentially 
 5130  be handled by the techniques described at the 
 5138  beginning of this section; but some of these 
 5146  distributions are so important in practice, they 
 5153  deserve special mention here.|'{A12}|\|ε(|↔O)|9|\The 
 5158  geometric distribution.|9|4|πIf some event occurs 
 5163  with probability |εp, |πthe number |εN |πof independent 
 5171  trials needed until the event _rst occurs (or 
 5179  between occurrences of the event) has the geometric 
 5187  distribution. We have |εN|4α=↓|41 |πwith probability 
 5193  |εp, N|4α=↓|42 |πwith probability (|ε1|4α_↓|4p)p,|4.|4.|4.|4
 5197  , N|4α=↓|4n |πwith probability (|ε1|4α_↓|4p)|gn|gα_↓|g1p. 
 5202  (|πThis is essentially the same situation as 
 5209  we have already considered in the gap test of 
 5218  Section 3.3.2; it is also directly related to 
 5226  the number of times certain loops in the algorithms 
 5235  of this section are executed, e.g., steps P1<P3 
 5243  of the polar method.)|'!|9|7A convenient way 
 5250  to generate a variable with this distribution 
 5257  is to set|'{A9}|εN|5|¬L|5|"p|πln|4|εU/|πln(|ε1|4α_↓|4p)|"P.|
 5260  E|;(34)|?{A9}|πTo check this formula, we observe 
 5268  that |"pln|4|εU/|πln(|ε1|4α_↓|4p)|"P|4α=↓|4n 
 5270  |πif and only if |εn|4α_↓|41|4|¬W|4|πln|4|εU/|πln(1|4α_↓|4|ε
 5274  p)|4|¬E|4n, |πthat is, (|ε1|4α_↓|4p)|gn|gα_↓|g1|4|¬Q|4U|4|¬R
 5277  |4(1|4α_↓|4p)|gn, |πand this happens with probability 
 5283  |εp(1|4α_↓|4p)|gn|gα_↓|g1 |πas required. Note 
 5287  thatln|4|εU |πcan be replaced by |→α_↓|εY, |πwhere 
 5294  |εY |πhas the exponential distribution with mean 
 5301  1.|'!|9|7The special case |εp|4α=↓|4|f1|d32|) 
 5306  can be handled more easily on a binary computer, 
 5315  since formula (34) becomes |εN|4|¬L|4|"p|→α_↓|πlog|β2|4|εU|¬
 5319  P, |πthat is, |εN |πis one more than the number 
 5329  of leading zero bits in the binary representation 
 5337  of {U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addiso
folio 173 galley 8
 5339  n-Wesley)|9f. 173|9ch. 3|9g. 8c|'{A6}|ε|*/(|↔P)|9|\The 
 5344  binomial distribution (t,|4p).|9|4|πIf some event 
 5349  occurs with probability |εP |πand if we carry 
 5357  out |εt |πindependent trials, the total number 
 5364  |εN |πof occurrences equals |εn |πwith probability 
 5371  (|ε|=|βn|gt)p|gn(1|4α_↓|4p)|gt|gα_↓|gn. (|πSee 
 5373  Section 1.2.10.) |πIn other words if we generate 
 5381  |εU|β1,|4.|4.|4.|4, U|βt, |πwe want to count 
 5387  how many of these are |→|¬W|εp. |πFor small |εt 
 5396  |πwe can obtain |εN |πin exactly this way.|'!|9|7For 
 5405  large |εt, |πwe can generate a beta variate |εX 
 5414  |πwith integer parameters |εa |πand |εb |πwhere 
 5421  |εa|4α+↓|4b|4α_↓|41|4α=↓|4t; |πthis e=ectively 
 5424  gives us the |εb|πth largest of |εt |πelements, 
 5432  without bothering to generate the other elements. 
 5439  Now if |εX|4|¬R|4p, |πwe set |εN|4|¬L|4N|β1 |πwhere 
 5446  |εN|β1 |πhas the binomial distribution (|εa|4α_↓|41,|4p/X), 
 5452  |πsince this tells us how many of |εa|4α_↓|41 
 5460  |πrandom numbers in the range [0,|4|εX) |πare 
 5467  |→|¬W|εp; |πand if |εX|4|¬W|4p, |πwe set |εN|4|¬L|4a|4α+↓|4N
 5473  |β1 |πwhere |εN|β1 |πhas the binomial distribution 
 5480  {H12}({H10}|εb|4α_↓|41,|4(p|4α_↓|4X)/(1|4α_↓|4X){H12}){H10},
 5480   |πsince |εN|β1 |πtells us how many of |εb|4α_↓|41 
 5489  |πrandom numbers in the range [|εX,|41) |πare 
 5496  |ε|→|¬Wp. |πBy choosing |εa|4α=↓|41|4α+↓|4|"lT/2|"L, 
 5500  |πthe parameter |εt |πwill be reduced to a reasonable 
 5509  size after about lg|4|εt |πreductions of this 
 5516  kind. (This approach is due to J. H. Ahrens, 
 5525  who has also suggested an alternative procedure 
 5532  for medium-sized |εt; |πsee exercise 27.)|'{A12}|ε|*/(|↔L)|9|
 5538  \The Poisson distribution |πwith mean |ε|≤m.|9|4|πThis 
 5544  distribution is related to the exponential distribution 
 5551  as the binomial distribution is related to the 
 5559  geometric: it represents the number of occurrences, 
 5566  per unit time, of an event which can occur at 
 5576  any instant of time; for example, the number 
 5584  of alpha particles emitted by radioactive substance 
 5591  in a single second has a Poisson distribution.|'
 5599  !|9|7|πAccording to this principle, we can produce 
 5606  a Poisson deviate |εN |πby generating independent 
 5613  exponential deviates |εX|β1,|4X|β2,|4.|4.|4. 
 5616  |πwith mean 1/|ε|≤m, |πstopping as soon as |εX|β1|4α+↓|4|¬O|
 5623  4|¬O|4|¬O|4α+↓|4X|βm|4|¬R|41; |πthen |εN|4|¬L|4m|4α_↓|41. 
 5626  |πThe probability that |εX|β1|4α+↓|4|¬O|4|¬O|4|¬O|4X|βm|4|¬R
 5629  |41 |πis the probability that a gamma deviate 
 5637  of order |εm |πis |ε|→|¬R|≤m, |πnamely |ε|↔j|ur|¬X|)|≤m|)|4t
 5643  |gm|gα_↓|g1e|gα_↓|gt|4dt/(m|4α_↓|41)*3; |πhence 
 5645  the probability that |εN|4α=↓|4n |πis|'{A9}|ε|(1|d2n*3|)|4|↔j
 5650  |ur|¬X|)|≤m|)|4t|gne|gα_↓|gt|4dt|4α_↓|4|(1|d2(n|4α_↓|41)*3|)|
 5650  4|↔j|ur|¬X|)|≤m|)|4t|gn|gα_↓|g1e|gα_↓|gt|4dt|4α=↓|4e|gα_↓|g|
 5650  ≤m|4|(|≤m|gn|d2n*3|)|4,!!n|4|¬R|40.|E|'(35)|?{A9}|πIf 
 5653  we generate exponential deviates by the logarithm 
 5660  method, the above recipe tells us to stop when 
 5669  |→α_↓(ln|4|εU|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4|πln|4|εU|βm)/|≤m
 5669  |4|¬R|41. |πSimplifying this expression, we see 
 5675  that the desired Poisson deviate can be obtained 
 5683  by calculating |εe|gα_↓|g|≤m, |πconverting it 
 5688  to a _xed-point representation, then generating 
 5694  one or more uniform deviates |εU|β1,|4U|β2,|4.|4.|4. 
 5700  |πuntil the product |εU|β1|4.|4.|4.|4U|βm|4|¬E|4e|gα_↓|g|≤m,
 5703   |πand _nally setting |εN|4|¬L|4m|4α_↓|41. |πOn 
 5709  the average this requires the generation of |ε|≤m|4α+↓|41 
 5717  |πuniform deviates, so it is a very useful approach 
 5726  when |ε|≤m |πis not too large. For very small 
 5735  |ε|≤m |πit would be better to generate exponential 
 5743  variables directly with a streamlined version 
 5749  of Algorithm S.|'!|9|7When |ε|≤m |πis large, 
 5756  we can obtain a method of order log|4|ε|≤m |πby 
 5765  using the fact that we know how to handle the 
 5775  gamma and binomial distributions for large orders: 
 5782  First generate |εX |πwith the gamma distribution 
 5789  of order |εm|4α=↓|4|"l|≤a|≤m|"L, |πwhere |ε|≤a 
 5794  |πis a suitable constant. (Since |εX |πis equivalent 
 5802  to |→α_↓ln(|εU|β1|4.|4.|4.|4U|βm), |πwe are essentially 
 5807  bypassing |εm |πsteps of the previous method.) 
 5814  If |εX|4|¬W|4|≤m, |πset |εN|4|¬L|4m|4α+↓|4N|β1, 
 5818  |πwhere |εN|β1 |πis a Poisson deviate of mean 
 5826  |ε|≤m|4α_↓|4X; |πand if |εX|4|¬R|4|≤m, |πset 
 5831  |εN|4|¬L|4N|β1, |πwhere |εN|β1 |πhas the binomial 
 5837  distribution (|εm|4α_↓|41,|4|≤m/X). |πThis method 
 5841  is due to J. H. Ahrens and U. Dieter, whose experiments 
 5852  suggest that |f7|d38|) is a good choice for |ε|≤a.|'
 5861  !|9|7|πThe validity of the above reduction when 
 5868  |εX|4|¬R|4|≤m |πis a consequence of the following 
 5875  important principle: ``Let |εX|β1,|4.|4.|4.|4, 
 5879  X|βm |πbe independent exponential deviates with 
 5885  the same mean; let |εS|βj|4α=↓|4X|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+
 5889  ↓|4X|βj |πand let |εV|βj|4α=↓|4S|βj/S|βm |πfor 
 5894  |ε1|4|¬E|4j|4|¬E|4m. |πThen the distribution 
 5898  of |εV|β1,|4V|β2,|4.|4.|4.|4, V|βm|βα_↓|β1 |πis 
 5902  the same as the distribution of |εm|4α_↓|41 |πindependent 
 5910  uniform deviates sorted into increasing order.'' 
 5916  To establish this principle formally, we compute 
 5923  the probability that |εV|β1|4|¬E|4v|β1,|4.|4.|4.|4, 
 5927  V|βm|βα_↓|β1|4|¬E|4v|βm|βα_↓|β1, |πgiven the 
 5930  value of |εS|βm|4α=↓|4s, |πfor arbitrarily values 
 5936  |ε0|4|¬E|4v|β1|4|¬E|4|¬O|4|¬O|4|¬O|4|¬Ev|βm|βα_↓|β1|4|¬E|41:
 5936  |'{A9}{M36}|↔j|urv|β1s|)0|)|4|≤m|πe|ε|gα_↓|gt|r1|g/|g|≤m|4dt
 5937  |β1|4|↔j|urv|β2sα_↓t|β1|)0|)|3|≤me|gα_↓|gt|r2|g/|g|≤m|4dt|β2
 5937  |4.|4.|4.|4|↔j|urv|βm|βα_↓|β1sα_↓t|β1|4α_↓|4|¬O|4|¬O|4|¬O|4α
 5937  _↓|4t|βm|βα_↓|β2|3|≤me|gα_↓|gt|rm|rα_↓|r1|g/|g|≤m|4dt|βm|βα_
 5937  ↓|β1|4|¬O|4|≤me|gα_↓|g(|gs|gα_↓|gt|r1|gα_↓|4|g|¬O|4|g|¬O|4|g
 5937  |¬O|4|gα_↓|gt|rm|rα_↓|r1|g)|g/|g|≤m|d2|↔j|urs|)0|)|3|ε|≤me|g
 5937  α_↓|gt|r1|g/|g|≤m|4dt|β1|4|↔j|ursα_↓t|β1|)0|)|3|≤me|gα_↓|gt|
 5937  rz|g/|g|≤m|4dt|β2|4.|4.|4.|3|↔j|ursα_↓t|β1α_↓|4|¬O|4|¬O|4|¬O
 5937  |4α_↓t|βm|βα_↓|β2|)0|)|3|≤me|gα_↓|gt|rm|rα_↓|r1|g/|g|≤m|4dt|
 5937  βm|βα_↓|β1|4|¬O|4|≤me|gα_↓|g(|gs|gα_↓|gt|r1|gα_↓|g|4|g|¬O|g|
 5937  4|g|¬O|g|4|g|¬O|g|4|gα_↓|gt|rm|rα_↓|r1|g)|g/|g|≤m|)|'
 5938  {A6}{M29}α=↓|4|(|↔j|urv|β1|)0|)|3du|β1|3|↔j|urv|β2|)u|β1|)|3
 5938  du|β2|4.|4.|4.|3|↔j|urv|βm|βα_↓|β1|)u|βm|βα_↓|β2|)|3du|βm|βα
 5938  _↓|β1|d2|↔j|ur1|)0|)|3|↔j|ur1|)u|β1|)|3du|β2|4.|4.|4.|3|↔j|u
 5938  r1|)u|βm|βα_↓|β2|)|3du|βm|βα_↓|β1|)|4,|?{A9}|πby 
 5940  making the substitution |εt|β1|4α=↓|4su|β1, t|β1|4α+↓|4t|β2|
 5944  4α=↓|4su|β2,|4.|4.|4.|4,|4t|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4t|β
 5944  m|βα_↓|β1|4α=↓|4su|βm|βα_↓|β1. |πThe latter ratio 
 5948  is the corresponding probability that uniform 
 5954  deviates |εU|β1,|4.|4.|4.|4, U|βm|βα_↓|β1 |πsatisfy 
 5958  |εU|β1|4|¬E|4v|β1,|4.|4.|4.|4, U|βm|βα_↓|β1|4|¬E|4v|βm|βα_↓|
 5959  β1, |πgiven that |εU|β1|4|¬E|4|¬O|4|¬O|4|¬O|4|¬E|4U|βm|βα_↓|
 5962  β1.|'!|9|7|πA more e∃cient but somewhat more 
 5969  complex technique for binomial and Poisson deviates 
 5976  is sketched in exercise 22.|'{A12}|≡G|≡.|9|≡F|≡o|≡r 
 5982  |≡f|≡u|≡r|≡t|≡h|≡e|≡r |≡r|≡e|≡a|≡d|≡i|≡n|≡g|≡.|9|4The 
 5984  book |εNon-Uniform Random Numbers |πby J. H. 
 5991  Ahrens and U. Dieter (New York: Wiley 1976) discusses 
 6000  many more algorithms for the generation of random 
 6008  variables with nonuniform distributions, together 
 6013  with a careful consideration of the e∃ciency 
 6020  of each technique on typical computers.|'!|9|7Exercise 
 6027  16 is recommended as a review of many of the 
 6037  techniques in this section.|'|Ha*?*?{U0}{H10L12M29}W58320#Comp
folio 176 galley 9
 6041  uter Programming|9(Knuth/Addison-Wesley)|9f. 
 6043  176|9ch. 3|9g. 9c|'{A6}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'
 6047  {A12}{H9L11}|5|5|≡1|≡.|9|4[|ε|*/|↔O|↔c|\] |πIf 
 6049  |ε|≤a |πand |ε|≤b |πare real numbers with |ε|≤a|4|¬W|4|≤b, 
 6057  |πhow would you generate a random real number 
 6065  uniformly distributed between |ε|≤a |πand |ε|≤b?|'
 6071  {A3}|5|5|≡2|≡.|9|4[|*/M16|\] |πAssuming that |εmU 
 6075  |πis a random integer between 0 and |εm|4α_↓|41, 
 6083  |πwhat is the |εexact |πprobability that |ε|"lkU|"L|4α=↓|4r,
 6089   |πif 0|4|¬E|4|εr|4|¬W|4k? |πCompare this with 
 6095  the desired probability |ε1/k.|'{A3}|5|5|≡3|≡.|9|4|ε[|*/|↔O|↔
 6099  M|\] |πDiscuss treating |εU |πas an integer and 
 6107  |εdividing |πby |εk |πto get a random integer 
 6115  between 0 and |εk|4α_↓|41, |πinstead of multiplying 
 6122  as suggested in the text. Thus (1) would be changed 
 6132  to|'{A9}|¬e|¬n|¬t|¬a!|¬o|;|¬l|¬d|¬x|6|6|¬u|;|¬d|¬i|¬v|6|6|¬k
 6135  |;{A6}with the result appearing in register X. 
 6143  Is this a good method?|'{A3}|5|5|≡4|≡.|9|4[|ε|*/M|↔P|↔c|\] 
 6149  |πProve the two relations in (7).|'{A3}|5|5|≡5|≡.|9|4[|ε|*/21
 6155  |\] |πSuggest an e∃cient method to compute a 
 6163  random variable with the distribution |εpx|4α+↓|4qx|g2|4α+↓|
 6168  4rx|g3, |πwhere |εp|4|¬R|40, q|4|¬R|40, r|4|¬R|40, 
 6173  |πand |εp|4α+↓|4q|4α+↓|4r|4α=↓|41.|'{A3}|5|56.|9|4[|ε|*/HM|↔P
 6175  |↔O|\] |πA quantity |εX |πis computed by the 
 6183  following method:|'{A9}!!|2``|εStep 1.|9|4|πGenerate 
 6187  two independent uniform deviates |εU, V.|'{A2}!!|9|3Step 
 6194  2.|9|4|πIf |εU|g2|4α+↓|4V|g2|4|¬R|41, |πreturn 
 6197  to step 1; otherwise set |εX|4|¬L|4U.''|'{A9}|πWhat 
 6204  is the distribution function of |εX? |πHow many 
 6212  times will step 1 be performed? (Give the mean 
 6221  and standard deviation.)|'{A3}|5|5|≡7|≡.|9|4[|ε|*/M|↔O|↔l|\] 
 6225  |πExplain why, in Marsaglia's method for normal 
 6232  deviates, the decision to make |εp|βj |πa multiple 
 6240  of |f1|d3256|) implies that we should take |εp|βj|4α=↓|4|"l6
 6247  4f|1(j/4)|"L/256, 1|4|¬E|4j|4|¬E|412.|'{A3}|5|5|≡8|≡.|9|4[|*/
 6249  |↔O|↔c|\] |πWhy are there ``skinny'' rectangles 
 6255  |εf|β1|β3,|4.|4.|4.|4,|4f|β2|β4 |πin Marsaglia's 
 6258  method, as well as the larger rectangles |εf|β1,|4.|4.|4.|4,
 6265  |4f|β1|β2? (|πWhy is this better than putting 
 6272  (|εf|β1,|4f|β1|β3), (f|β2,|4f|β1|β4),|4.|4.|4. 
 6274  |πtogether into single large rectangles?)|'{A3}|5|5|≡9|≡.|9|
 6279  4[|ε|*/HM|↔O|↔c|\] |πWhy is the curve |εf|1(x) 
 6285  |πof Fig. 9 concave downward for |εx|4|¬W|41, 
 6292  |πconcave upward for |εx|4|¬Q|41?|'{A3}|≡1|≡0|≡.|9|4[|ε|*/HM|
 6296  ↔P|↔L|\] |πGive formulas for the numbers |εa|βj, 
 6303  b|βj, p|βj|βα+↓|β2|β4 |πsuch that the densities 
 6309  |εf|βj|βα+↓|β2|β4(x) |πof Fig. 9 can be calculated 
 6316  by Algorithm L with |εa|4α=↓|4a|βj, b|4α=↓|4b|βj, 
 6322  h|4α=↓|4|f1|d34|), |πand |εs|4α=↓|4|f1|d34|)(j|4α_↓|41), 
 6325  |πfor |ε1|4|¬E|4j|4|¬E|412.|'{A3}|≡1|≡1|≡.|9|4[|*/HM|↔P|↔p|\]
 6327   |πProve that steps M8<M9 of Algorithm M generate 
 6336  a random variable with the appropriate tail of 
 6344  the normal distribution; i.e., the probability 
 6350  that |εX|4|¬E|4x |πshould be|'{A9}|ε|↔j|urx|)3|)|3e|gα_↓|gt|
 6354  i2|g/|g2|4dt|↔h|2|↔j|ur|¬X|)3|)|3e|gα_↓|gt|i2|g/|g2|4dt,!!x|
 6354  4|¬R|43.|;{A9}|πwhen |εx|4|¬R|43.|'{A3}|≡1|≡2|≡.|9|4[|*/HM|↔M
 6357  |↔o|\] |πTeichroew's polynomial, Eq. (22), is 
 6363  a truncated sum of Chebyshev polynomials and 
 6370  it may not be the best possible polynomial approximation 
 6379  of degree 9|'{A12}[|εHint|*/:|9|\|πIt is a special 
 6386  case of the rejection method, with |εg(t)|4α=↓|4Cte|gα_↓|gt|
 6392  i2|g/|g2 |πfor some |εC.]|'{A3}|≡1|≡2|≡.|9|4[|ε|*/HM|↔P|↔L|\]
 6396   |π(R. P. Brent.) Prove that the numbers |εa|βj 
 6405  |πde_ned in (23) satisfy |εa|=|βj|g2/2|4α_↓|4a|ur2|)jα_↓1|)/
 6409  2|4|¬W|4|πln|42 for all |εj|4|¬R|41. [Hint|*/: 
 6414  |\|πIf |εf|1(x)|4α=↓|4e|gx|i2|g/|g2|4|¬J|ur|¬X|)x|)|4e|gα_↓|
 6415  gt|i2|g/|g2|4dt, |πshow that |εf|1(x)|4|¬W|4f|1(y) 
 6419  |πfor |ε0|4|¬E|4x|4|¬W|4y.]|'{A3}|≡1|≡3|≡.|9|4[|ε|*/HM|↔P|↔C|
 6421  \] |πGiven a set of |εn |πindependent normal 
 6429  deviates, |εX|β1, X|β2,|4.|4.|4.|4,|4X|βn, |πwith 
 6433  mean 0 and variance 1, show how to _nd constants 
 6443  |εb|βj |πand |εa|βi|βj, 1|4|¬E|4j|4|¬E|4i|4|¬E|4n, 
 6447  |πso that if|'{A9}!!|εY|β1|4α=↓|4b|β1|4α+↓|4a|β1|β1X|β1,!!Y|
 6450  β2|4α=↓|4b|β2|4α+↓|4a|β2|β1X|β1|4α+↓|4a|β2|β2X|β2,!!.|4.|4.|
 6450  4,|'{A4}|εY|βn|4α=↓|4b|βn|4α+↓|4a|βn|β1X|β1|4α+↓|4a|βn|β2X|β
 6451  2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓a|βn|βnX|βn,!!|?{A9}|πthen 
 6453  |εY|β1, Y|β2,|4.|4.|4.|4,|4Y|βn |πare dependent, 
 6457  normally distributed variables, |εY|βj |πhas 
 6462  mean |ε|≤m|βj, |πand the |εY'|πs have a given 
 6470  covariance matrix (|εc|βi|βj). (|πThe covariance, 
 6475  |εc|βi|βj, |πof |εY|βi |πand |εY|βj |πis de_ned 
 6482  to be the average value of (|εY|βi|4α_↓|4|≤m|βi)(Y|βj|4α_↓|4
 6488  |≤m|βj). |πIn particular, |εc|βj|βj |πis the 
 6494  variance of |εY|βj, |πthe square of its standard 
 6502  deviation. Not all matrices (|εc|βi|βj) |πcan 
 6508  be covariance matrices, and your construction 
 6514  is, of course, only supposed to work whenever 
 6522  a solution to the given conditions is possible.)|'
 6530  {A3}|≡1|≡4|≡.|9|4[|ε|*/M|↔P|↔O|\] |πIf |εX |πis 
 6534  a random variable with continuous distribution 
 6540  |εF(x), |πand if |εc |πis a constant, what is 
 6549  the distribution of |εcX?|'{A3}|≡1|≡5|≡.|9|4[|*/HM21|\] 
 6554  |πIf |εX|β1 |πand |εX|β2 |πare independent random 
 6561  variables with the respective distributions |εF|β1(x) 
 6567  |πand |εF|β2(x), |πand the densities |εf|β1(x)|4α=↓|4F|=|β1|
 6572  ¬S(x), f|β2(x)|4α=↓|4F|=|β2|¬S(x), |πwhat are 
 6576  the distribution and density functions of the 
 6583  quantity |εX|β1|4α+↓|4X|β2?|'{A3}|≡1|≡6|≡.|9|4[|ε|*/HM22|\] 
 6586  |π(J. H. Ahrens.) Develop an algorithm for gamma 
 6594  deviates of order |εa |πwhen |ε0|4|¬W|4|εa|4|¬E|41, 
 6600  |πusing the rejection method with |εcg(t)|4α=↓|4t|ga|gα_↓|g1
 6605  /|≤)(a) |πfor |ε0|4|¬W|4t|4|¬W|41, |εe|gα_↓|gt/|≤)(a) 
 6609  |πfor |εt|4|¬R|41.|'{A3}|≡1|≡7|≡.|9|4[|ε|*/M|↔P|↔M|\] 
 6612  |πWhat is the |εdistribution function F(x) |πfor 
 6619  the geometric distribution with probability |εp? 
 6625  |πWhat is the |εgenerating function G(z)? |πWhat 
 6632  are the mean and standard deviation of this distribution?|'
 6641  {A3}|≡1|≡8|≡.!|9|7[|ε|*/N|↔P|↔M|\] |πSuggest a 
 6644  method to compute a random integer |εN |πfor 
 6652  which |εN |πtakes the value |εn |πwith probability 
 6660  |εnp|g2(1|4α_↓|4p)|gn|gα_↓|g1, n|4|¬R|40. (|πThe 
 6663  case of particular interest is when |εp |πis 
 6671  rather small.)|'{A3}|≡1|≡9|≡.|9|4[|ε|*/|↔P|↔P|\] 
 6674  |πThe |εnegative binomial distribution (t,|4p) 
 6679  |πhas integer values |εN|4α=↓|4n |πwith probability 
 6685  (|ε|ft|2α_↓|21|2α+↓|2n|d5n|))p|gt(1|4α_↓|4p)|gn. 
 6686  (|πUnlike the ordinary binomial distribution, 
 6691  |εt |πneed not be an integer, since this quantity 
 6700  is nonnegative for all |εn |πwhenever |εt|4|¬Q|40.) 
 6707  |πGeneralizing exercise 18, explain how to generate 
 6714  integers |εN |πwith this distribution when |εt 
 6721  |πis a small positive integer. What method would 
 6729  you suggest if |εt|4α=↓|4p|4α=↓|4|f1|d32|)?|'
 6733  {A3}|≡2|≡0|≡.|9|4[|ε|*/|↔P|↔L|\] |π(``Marsaglia 
 6735  tables.'') Suppose we want to compute a random 
 6743  variable |εX |πwith the following _nite distribution:|'
 6750  {A9}|hWith|6probability|4|∂α=↓|4|f1|d3512|)!|f8|d3512|)!|f3|
 6750  d3512|)!|f8|d3512|)!|f3|d3512|)!|f6|d3512|)|E|n|;
 6751  | |πValue|6of|6|εX|4|Lα=↓|4|5x|β1|6!|5x|β2|6!|5x|β3|6!|5x|β4
 6751  |6!|5x|β5|6!|5x|β6>{A4}| |πWith|6probability|4|Lα=↓|4|f90|d3
 6752  512|)!|f81|d3512|)!|f131|d3512|)!|f10|d3512|)!|f32|d3512|)!|
 6752  f168|d3512|)>{A9}|πShow how the value of |εX 
 6759  |πcan be e∃ciently obtained from a random nine-bit 
 6767  binary integer, |εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8b|β9, 
 6770  |πusing a method analogous to (16) and using 
 6778  auxiliary tables containing less than 35 entries 
 6785  in all.|'|Ha*?*?*?*?{U0}{H10L12M29}W58320#Computer 
folio 178 galley 10
 6788  Programming|9(Knuth/Addison-Wesley)|9f. 178|9ch. 
 6790  3|9g. 10c|'{A6}{H9L11}|π|≡2|≡1|≡.|9|4[|ε|*/|↔P|↔M|\] 
 6793  |πThe situation in the preceding exercise is 
 6800  idealized in the sense that we usually will not 
 6809  have all probabilities an exact multiple of a 
 6817  convenient number like |f1|d3512|). Suppose the 
 6823  probabilities which appear in that exercise are 
 6830  modi_ed so that they are actually equal to the 
 6839  following:|'{A9}{A16}{M32}With|4probability|4|∂α=↓|4|f89|d35
 6840  12|4α+↓|4|≤e|β1!|∂|f80|d3512|)|4α+↓|4|≤e|β2!|∂|f131|d3512|)|
 6840  4α+↓|4|≤e|β3!|∂|f10|d3512|)|4α+↓|4|≤e|β4!|∂|f32|d3512|)|4α+↓
 6840  |4|≤e|β5!|∂|f168|d3512|)|4α+↓|4|≤e|β6|E|'{B16}| Value|4of|4|
 6841  εX|4|Lα=↓!|9x|β1|L!|9x|β2|L!|9x|β3|L!|9x|β4|L!|9x|β5|L!|9x|β
 6841  6>{A25}{M29}|πHere |ε0|4|¬E|4|≤e|βj|4|¬W|4|f1|d3512|), 
 6844  |πand |ε|≤e|β1|4α+↓|4|≤e|β2|4α+↓|4|≤e|β3|4α+↓|4|≤e|β4|4α+↓|4
 6845  |≤e|β5|4α+↓|4|≤e|β6|4α=↓|4|f2|d3512|). |πShow 
 6847  how to make this choice e∃ciently, when given 
 6855  a uniform binary deviate |εU|4α=↓|4.|εb|β1b|β2|4.|4.|4.|4b|β
 6859  m |πwhere |εm |πis relatively large. As in exercise 
 6868  20, try to do this using auxiliary tables containing 
 6877  less than 35 entries in all.|'{A3}|≡2|≡2|≡.|9|4[|ε|*/HM|↔P|↔P
 6883  |\] |πCan the exact Poisson distribution for 
 6890  large |ε|≤m |πbe obtained by generating an appropriate 
 6898  normal deviate, converting it to an integer in 
 6906  some convenient way, and applying a (possibly 
 6913  complicated) correction a small percent of the 
 6920  time?|'{A3}|≡2|≡3|≡.|9|4[|ε|*/HM|↔P|↔L|\] (|πJ. 
 6923  von Neumann.) Are the following two ways to generate 
 6932  a random quantity |εX |πequivalent (i.e., does 
 6939  the quantity |εX |πhave the same distribution)?|'
 6946  {A3}!!|2|εMethod|6|ε|*/|↔O|\:|5|5|πSet |εX|5|¬L|5|πsin{H12}({
 6947  H10}(|ε|≤p/2)U{H12}){H10}, |πwhere |εU |πis uniform.|'
 6952  {A3}{I1.7l}|εMethod|6|*/|↔P|\:|5|5|πGenerate two 
 6954  uniform deviates, |εU, V, |πand if |εU|g2|4α+↓|4V|g2|4|¬R|41
 6960  , |πrepeat until |εU|g2|4α+↓|4V|g2|4|¬W|41. |πThen 
 6965  set |εX|5|¬L|5|¬GU|g2|4α_↓|4V|g2|¬G/(U|g2|4α+↓|4V|β2).|'
 6967  {A3}{IC}|π|≡2|≡4|≡.|9|4[|ε|*/HM|↔M|↔c|\] (|πS. 
 6969  Ulam, J. von Neumann.) Let |εV|β0 |πbe a randomly 
 6978  selected real number between 0 and 1, and de_ne 
 6987  the sequence |¬D|εV|βn|¬F |πby the rule |εV|βn|βα+↓|β1|4α=↓|
 6993  44V|βn(1|4α_↓|4V|βn). |πNow if this computation 
 6998  is done with perfect accuracy, the result should 
 7006  be a sequence with the distribution sin|g2|4|ε|≤pU, 
 7013  |πwhere |εU |πis uniform, i.e., with distribution 
 7020  function|'{A9}|εF(x)|4α=↓|3|↔j|urx|)0|)|3dx/{H11}|¬H{H10}|v2
 7021  2|≤px(1|4α_↓|4x)|).|;{A9}|πFor if we write |εV|βn|4α=↓|4|πsi
 7026  n|g2|4|ε|≤pU|βn, |πwe _nd that |εU|βn|βα+↓|β1|4α=↓|4(2U|βn)|
 7030  πmod 1; and by the fact that almost all real 
 7040  numbers have a random binary expansion (see Section 
 7048  3.5), this sequence |εU|βn |πis equidistributed. 
 7054  But if the computation of |εV|βn |πis done with 
 7063  only _nite accuracy, the above argument breaks 
 7070  down because we soon are dealing with noise from 
 7079  the roundo= error. [|εReference|*/: |\|πvon Neumann's 
 7085  |εCollected Works |π|≡5|≡, 768<770.]|'!|9|7Analyze 
 7090  the sequence |¬D|εV|βn|¬F |πde_ned above when 
 7096  only _nite accuracy is present, both empirically 
 7103  (for various di=erent choices of |εV|β0) |πand 
 7110  theoretically. Does the sequence have a distribution 
 7117  resembling the expected distribution?|'{A3}|≡2|≡5|≡.|9|4[|ε|
 7121  */M|↔P|↔C|\] |πLet |εX|β1, X|β2,|4.|4.|4.|4,|4X|β5 
 7125  |πbe binary words each of whose bits is independently 
 7134  0 or 1 with probability |f1|d32|). What is the 
 7143  probability that a given bit position of |εX|β1|4|↔I|4{H12}(
 7150  {H10}X|β2|4|↔i|4{H12}({H10}X|β3|4|↔I|4(X|β4|4|↔i|4X|β5){H12}
 7150  )){H10} |πcontains a 1? Generalize.|'{U31}1|'
 7156  {H9L11M29}{A3}|π|≡2|≡6|≡.|9|4[|ε|*/M|↔O|↔l|\] 
 7157  |πLet |εN|β1 |πand |εN|β2 |πbe independent Poisson 
 7164  deviates with respective means |ε|≤m|β1 |πand 
 7170  |ε|≤m|β2, |πwhere |ε|≤m|β1|4|¬Q|4|≤m|β2|4|¬R|40. 
 7173  |πProve or disprove: (a) |εN|β1|4α+↓|4N|β2 |πhas 
 7179  the Poisson distribution with mean |ε|≤m|β1|4α+↓|4|≤m|β2. 
 7185  |π(b) |εN|β1|4α_↓|4N|β2 |πhas the Poisson distribution 
 7191  with mean |ε|≤m|β1|4α_↓|4|≤m|β2.|'{A3}|≡2|≡7|≡.|9|4[|ε|*/|↔P|
 7194  ↔P|\] (|πJ. H. Ahrens.) On most binary computers 
 7202  there is an e∃cient way to count the number of 
 7212  1's in a binary word (cf. Section 7.1). Hence 
 7221  there is a nice way to obtain the binomial distribution 
 7231  (|εt,|4p) |πwhen |εp|4α=↓|4|f1|d32|), |πsimply 
 7235  by generating |εt |πrandom bits and counting 
 7242  the number of |ε|*/|↔O|\|π's.|'!|9|7|πDesign an 
 7248  algorithm which produces the binomial distribution 
 7254  (|εt,|4p) |πfor |εarbitrary p, |πusing only a 
 7261  subroutine for the special case |εp|4α=↓|4|f1|d32|) 
 7267  |πas a source of random data. [|εHint: |πSimulate 
 7275  a process which _rst looks at the most signi_cant 
 7284  bits of |εt |πuniform deviates, then at the second 
 7293  bit of those deviates whose leading bit is not 
 7302  su∃cient to determine whether or not their value 
 7310  is |→|¬W|εp, |πetc.]|'{A3}|≡2|≡8|≡.|9|4[|ε|*/|↔P|↔C|\] 
 7314  (|πA. J. Walker.) The _nite distribution in exercise 
 7322  20 could also be generated by the following method: 
 7331  ``Generate a uniform deviate |εU, |πand set |εN|5|¬L|5|"l6U|
 7338  "L. |πIf |εU|4|¬W|4P|βN |πset |εX|4|¬L|4x|βN|βα+↓|β1, 
 7343  |πotherwise set |εX|4|¬L|4Q|βN, |πwhere (|εP|β0,|4.|4.|4.|4,
 7347  |4P|β5)|4α=↓|4(|f256|d31536|), |f499|d31536|), 
 7349  |f745|d31536|), |f798|d31536|), |f1120|d31536|), 
 7352  |f1535|d31536|)) |πand (|εQ|β1,|4.|4.|4.|4,|4Q|β5)|4α=↓|4(x|
 7354  β1,|4x|β6,|4x|β6,|4x|β3,|4x|β1).'' |πShow that 
 7357  the general _nite distribution (2) can always 
 7364  be handled in a similar way, with appropriate 
 7372  tables (|εP|β0,|4.|4.|4.|4,|4P|βk|βα_↓|β1), (Q|β0,|4.|4.|4.|
 7374  4,|4Q|βk|βα_↓|β1).|'{A3}|≡2|≡9|≡.|9|4[|ε|*/HM|↔L|↔C|\] 
 7376  (|πR. P. Brent.) Develop a method to generate 
 7384  a random point on the surface of the ellipsoid 
 7393  de_ned by |¬K|4a|βkx|ur2|)k|)|4α=↓|41, |πwhere 
 7397  |ε|εa|β1|4|¬R|4|¬O|4|¬O|4|¬O|4|¬R|4a|βn|4|¬Q|40.|'
 7398  |H{A6}{H10L12M29}|πThe (|εt|4α+↓|41)|πst record 
 7401  should be selected with probability (|εn|4α_↓|4m)/(N|4α_↓|4t
 7406  ), |πif |εm |πitems have already been selected. 
 7414  This is the appropriate probability, since of 
 7421  all the possible ways to choose |εn |πthings 
 7429  from |εN |πsuch that |εm |πvalues occur in the 
 7438  _rst |εt, |πexactly|'{A9}|ε|↔a|(N|4α_↓|4t|4α_↓|41|d5n|4α_↓|4
 7441  m|4α_↓|41|)|↔s|↔h|↔a|(N|4α_↓|4t|d5n|4α_↓|4m|)|↔s|4α=↓|4|(n|4
 7441  α_↓|4m|d2N|4α_↓|4t|)|E|;(1)|?{A9}|πof these select 
 7446  the (|εt|4α+↓|41)|πst element.|'!|9|7The idea 
 7451  developed in the preceding paragraph leads immediately 
 7458  to the following algorithm:|'|Hβ{U0}{H10L12M29}W58320#Comput
folio 182 galley 11
 7462  er Programming|9(Knuth/Addison-Wesley)|9f. 182|9ch. 
 7465  3|9g. 11c|'{A6}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡S 
 7469  (|εSelection sampling technique).|9|4|πTo select 
 7473  |εn |πrecords at random from a set of |εN, |πwhere 
 7483  |ε0|4|¬W|4n|4|¬E|4N.|'{A3}{I1.7H}|≡S|≡1|≡.|9|π[Initialize.] 
 7485  Set |εt|4|¬L|40, m|4|¬L|40. (|πDuring this algorithm, 
 7491  |εm |πrepresents the number of records selected 
 7498  so far, and |εt |πis the total number of input 
 7508  records we have dealt with.)|'{A3}|≡S|≡2|≡.|9[|πGenerate 
 7514  |εU.] |πGenerate a random number |εU, |πuniformly 
 7521  distributed between zero and one.|'{A3}|≡S|≡3|≡.|9[Test.] 
 7527  If (|εN|4α_↓|4t)U|4|¬R|4n|4α_↓|4m, |πgo to step 
 7532  S5.|'{A3}|≡S|≡4|≡.|9[Select.] Select the next 
 7537  record for the sample, and increase |εm |πand 
 7545  |εt |πby 1. If |εm|4|¬W|4n, |πgo to step S2; 
 7554  otherwise the sample is complete and the algorithm 
 7562  terminates.|'{A3}|≡S|≡5|≡.|9[|πSkip.] Skip the 
 7566  next record (do not include it in the sample), 
 7575  increase |εt |πby 1, and go to step S2.|'{A12}{IC}|π!|9|7Thi
 7584  s algorithm, at _rst glance, may appear to be 
 7593  unreliable and, in fact, incorrect; but a careful 
 7601  analysis (see the exercises below) shows that 
 7608  it is completely trustworthy. It is not di∃cult 
 7616  to verify that|'{A3}a)|5|5At most |εN |πrecords 
 7623  are input (we never run o= the end of the _le 
 7634  before choosing |εn |πitems).|'{A3}b)|9The sample 
 7640  is completely unbiased; in particular, the probability 
 7647  that any given element is selected, e.g., the 
 7655  last element of the _le, is |εn/N.|'{A3}!|9|7|πStatement 
 7663  (b) is true in spite of the fact that we are 
 7674  |εnot |πselecting the (|εt|4α+↓|41)|πst item 
 7679  with probability |εn/N, |πwe select it with the 
 7687  probability in Eq. (1)*3 This has caused some 
 7695  confusion in the published literature. Can the 
 7702  reader explain this seeming contradiction?|'!|9|7(|εNote|*/: 
 7708  |\|πWhen using Algorithm S, care should be taken 
 7716  to use a di=erent source of random numbers |εU 
 7725  |πeach time the program is run, to avoid connections 
 7734  between the samples obtained on di=erent days. 
 7741  This can be done, for example, by choosing a 
 7750  di=erent value of |εX|β0 |πfor the linear congruential 
 7758  method each time; |εX|β0 |πcould be set to the 
 7767  current date, or to the last |εX |πvalue generated 
 7776  on the previous run of the program.)|'!|9|7We 
 7784  will usually not have to pass over all |εN |πrecords; 
 7794  in fact, since (b) above says the alst record 
 7803  is selected with probability |εn/N, |πwe will 
 7810  terminate the algorithm |εbefore |πconsidering 
 7815  the last record exactly (1|4α_↓|4|εn/N) |πof 
 7821  the time. The average number of records considered 
 7829  when |εn|4α=↓|42 |πis about |f2|d33|)|εN, |πand 
 7835  the general formulas are given in exercises 5 
 7843  and 6.|'!|9|7A problem arises if we don't know 
 7852  the value of |εN |πin advance, since the precise 
 7861  value of |εN |πis crucial in Algorithm S. Suppose 
 7870  we want to select |εn |πitems at random from 
 7879  a _le, without knowing exactly how many are present 
 7888  in that _le. We could _rst go through and count 
 7898  the records, then take a second pass to select 
 7907  them; but it is generally better to sample |εm|4|¬R|4n 
 7916  |πof the original items on the _rst pass, where 
 7925  |εm |πis much less than |εN, |πso that only |εm 
 7935  |πitems must be considered on the second pass. 
 7943  The trick, of course, is to do this in such a 
 7954  way that the _nal result is a truly random sample 
 7964  of the original _le.|'!|9|7|πSince we don't know 
 7972  when the input is going to end, we must keep 
 7982  track of a random sample of the input records 
 7991  seen so far, thus always being prepared for the 
 8000  end. As we read the input we will construct a 
 8010  ``reservoir'' which contains only those |εm |πrecords 
 8017  which have appeared among the previous samples. 
 8024  The _rst |εn |πrecords always go into the reservoir. 
 8033  When the (|εt|4α+↓|41)|πst record is being input, 
 8040  for |εt|4|¬E|4n, |πwe will have in memory a table 
 8049  of |εn |πindices pointing to those records in 
 8057  the reservoir which belong to the random sample 
 8065  we have chosen from the _rst |εt |πrecords. The 
 8074  problem is to maintain this situation with |εt 
 8082  |πincreased by one, namely to _nd a new random 
 8091  sample from among the |εt|4α+↓|41 |πrecords now 
 8098  known to be present. It is not hard to see that 
 8109  we should include the new record in the new sample 
 8119  with probability |εn/(t|4α+↓|41), |πand in such 
 8125  a case it should replace a random element of 
 8134  the previous sample.|'!|9|7Thus, the following 
 8140  procedure does the job:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
 8145  |≡R (|εReservoir sampling).|9|4|πTo select |εn 
 8150  |πrecords at random from a _le of unknown size 
 8159  |→|¬R|εn, |πgiven |εn|4|¬Q|40. |πAn auxiliary 
 8164  _le called the ``reservoir'' contains all records 
 8171  which are condidates for the _nal sample. The 
 8179  algorithm uses a table of distinct indices |εI[j], 
 8187  1|4|¬E|4j|4|¬E|4n, |πeach of which points to 
 8193  one of the records in the reservoir.|'{A3}{I1.8H}|π|≡R|≡1|≡.
 8200  |9[Initialize.] Input the _rst |εn |πrecords 
 8206  and copy them to the reservoir. Set |εI[j]|4|¬L|4j 
 8214  |πfor |ε1|4|¬E|4j|4|¬E|4n, |πand set |εt|5|¬L|5m|5|¬L|5n. 
 8219  |πDuring this algorithm, |¬T|εI[1],|4.|4.|4.|4,|4I[n]|¬S 
 8223  |πpoint to the records in the current sample, 
 8231  |εm |πis the size of the reservoir, and |εt |πis 
 8241  the number of input records dealt with so far.)|'
 8250  {A3}|≡R|≡2|≡.|9[End of _le?] If there are no 
 8257  more records to be input, go to step R6.|'{A6}|≡R|≡3|≡.|9[Ge
 8266  nerate and test.] Increase |εt |πby 1, then generate 
 8275  a random integer |εM |πbetween 1 and |εt (|πinclusive). 
 8284  If |εM|4|¬Q|4n, |πgo to R5.|'{A3}|≡R|≡4|≡.|9[Add 
 8290  to reservoir.] Copy the next record of the input 
 8299  _le to the reservoir, increase |εm |πby 1, and 
 8308  set |εI[M]|5|¬L|5m. (|πThe record previously 
 8313  pointed to by |εI[M] |πis being replaced in the 
 8322  sample by the new record.) Go back to R2.|'{A3}|≡R|≡5|≡.|9[S
 8331  kip.] Skip over the next record of the input 
 8340  _le (do not include it in the reservoir), and 
 8349  return to step R2.|'{A3}|≡R|≡6|≡.|9[Second pass.] 
 8355  Sort the |εI |πtable entries so that |εI[1]|4|¬W|4|¬O|4|¬O|4
 8362  |¬O|4I[n]; |πthen go through the reservoir, copying 
 8369  the records with these indices into the output 
 8377  _le which is to hold the _nal sample.|'{A12}{IC}!|9|7Algorit
 8385  hm R is due to Alan G. Waterman. The reader may 
 8396  wish to work out the example of its operation 
 8405  which appears in exercise 9.|'!|9|7If the records 
 8413  are su∃ciently short, it is of course unnecessary 
 8421  to have a reservoir at all; we can keep the |εn 
 8432  |πrecords of the current sample in memory at 
 8440  all times (see exercise 10).|'!|9|7|πThe natural 
 8447  question to ask about Algorithm R is, ``What 
 8455  is the expected size, of the reservoir?'' Exercise 
 8463  11 shows that the average value of |εm |πis exactly 
 8473  |εn(1|4α+↓|4H|βN|4α_↓|4H|βn); |πthis is approximately 
 8477  |εn{H12}({H10}1|4α+↓|4|πln(|εN/n){H12}){H10}. 
 8478  So if |εN/n|4α=↓|41000, |πthe reservoir will 
 8484  contain only a{U0}{H10L12M29}W58320#Computer 
folio 186 galley 12
 8487  Programming|9(Knuth/Addison-Wesley)|9f. 186|9ch. 
 8489  3!g. 12c|'{A6}!|9|7|πNote that Algorithms S and 
 8496  R can be used to obtain samples for several independent 
 8506  categories simultaneously. For example, if we 
 8512  have a large _le of names and addresses of U.S. 
 8522  residents, we could pick random samples of exactly 
 8530  10 people from each of the 50 states without 
 8539  making 50 passes through the _le, and without 
 8547  _rst sorting the _le by state.|'!|9|7Algorithm 
 8554  S and a number of other sampling techniques are 
 8563  discussed in a paper by C. T. Fan, Mervin E. 
 8573  Muller, and Ivan Rezucha, |εJournal of the American 
 8581  Statistical Association |π|≡5|≡7 (1962), 387<402. 
 8586  The method was independently discovered by T. 
 8593  G. Jones, |εCACM |π|≡5 (1962), 343.|'!|9|7The 
 8600  |εsampling problem, |πas described here, can 
 8606  be regarded as the computation of a random |εcombination, 
 8615  |πaccording to the conventional de_nition of 
 8621  combinations of |εN |πthings taken |εn |πat a 
 8629  time (see Section 1.2.6). Now let us consider 
 8637  the problem of computing a random |εpermutation 
 8644  |πof |εt |πobjects; we will call this the |εshu+ing 
 8653  problem, |πsince shu+ing a deck of cards is nothing 
 8662  more than subjecting it to a random permutation.|'
 8670  !|9|7A moment's re⊗ection is enough to convince 
 8677  oneself that the approaches people traditionally 
 8683  use to shu+e cards are miserably inadequate; 
 8690  there is no hope of obtaining each of the |εt*3 
 8700  |πpermutations with anywhere near equal probability 
 8706  by such methods. It has been said that expert 
 8715  bridge players make use of this fact when deciding 
 8724  whether or not to ``_nesse.''|'!|9|7|πWhen |εt 
 8731  |πis small, we can obtain a random permutation 
 8739  very quickly by generating a random integer between 
 8747  1 and |εt*3. |πFor example, when |εt|4α=↓|44, 
 8754  |πa random number between 1 and 24 su∃ces to 
 8763  select a random permutation from a table of all 
 8772  possibilities. But for alrge |εt, |πit is necessary 
 8780  to be more careful if we want to claim that each 
 8791  permutation is equally likely, since |εt*3 |πis 
 8798  much larger than the accuracy of individual random 
 8806  numbers.|'!|9|7A suitable shu+ing procedure can 
 8812  be obtained by recalling Algorithm 3.3.2P, which 
 8819  gives a simple correspondence between each of 
 8826  the |εt*3 |πpossible permutations and a sequence 
 8833  of numbers (|εc|β1,|4c|β2,|4.|4.|4.|4,|4c|βt|βα_↓|β1) 
 8836  |πwith |ε0|4|¬E|4c|βj|4|¬E|4j. |πIt is easy to 
 8842  compute such a set of numbers at random, and 
 8851  we can use the correspondence to produce a random 
 8860  permutation.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
 8862  |≡P (|εShu∪ing).|9|4|πLet |εX|β1, X|β2, |4.|4.|4.|4,|4X|βt 
 8867  |πbe a set of |εt |πnumbers to be shu+ed.|'{A3}{I1.8H}|≡P|≡1
 8876  |≡.|9[Initialize.] Set |εj|5|¬L|5t.|'{A3}|π|≡P|≡2|≡.|9[Gener
 8879  ate |εU.] |πGenerate a random number U, |πuniformly 
 8887  distributed between zero and one.|'{A3}|≡P|≡3|≡.|9[Exchange.
 8892  ] Set |εk|5|¬L|5|"ljU|"L|4α+↓|41. (|πNow |εk 
 8897  |πis a random integer between 1 and |εj). |πExchange 
 8906  |εX|βk|5|≠m|5X|βj.|'{A3}|π|≡P|≡4|≡.|9[Decrease 
 8908  |εj.] |πDecrease |εj |πby 1. If |εj|4|¬Q|41, 
 8915  |πreturn to step P2.|'{A12}{IC}!|9|7This algorithm 
 8921  was _rst published by L. E. Moses and R. V. Oakford, 
 8932  in |εTables of Random Permutations (|πStanford 
 8938  University Press, 1963), and by R. Durstenfeld, 
 8945  |εCACM |π|≡7 (1964), 420. It can be used also 
 8954  to obtain random combinations (see exercise 14).|'
 8961  {A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'{A12}{H9L11}|5|5|≡1|≡.|9|4
 8962  [|ε|*/M|↔O|↔P|\] |πExplain Eq. (1).|'{A3}|5|5|≡2|≡.|9|4[|ε|*/|
 8966  ↔P|↔c|\] |πProve that Algorithm S never tries 
 8973  to read more than |εN |πrecords of its input 
 8982  _le.|'{A3}|5|5|≡3|≡.|9|4[|ε|*/|↔P|↔P|π|\] The 
 8985  (|εt|4α+↓|41)|πst item in Algorithm S is selected 
 8992  with probability (|εn|4α_↓|4m)/(N|4α_↓|4t), not 
 8996  n/N, |πyet the text claims that the sample is 
 9005  unbiased so each item should be selected with 
 9013  the |εsame |πprobability*3 How can both of these 
 9021  statements be true?|'{A3}|5|5|≡4|≡.|9|4[|ε|*/M|↔P|↔L|\] 
 9025  |πLet |εp(m,|4t) |πbe the probability that exactly 
 9032  |εm |πitems are selected from among the _rst 
 9040  |εt |πin the selection sampling technique. Show 
 9047  directly from Algorithm S that|'{A9}|εp(m,|4t)|4α=↓|4|↔a|(t|
 9052  d5m|)|↔s|↔a|(N|4α_↓|4t|d5n|4α_↓|4m|)|↔s|↔h|↔a|(N|d5n|)|↔s,!!
 9052  |πfor!!|ε0|4|¬E|4t|4|¬E|4N.|;{A9}|5|5|≡5|≡.|9|4[|ε|*/M|↔P|↔M|
 9053  \] |πWhat is the average value of |εt |πwhen 
 9062  Algorithm S terminates? (In other words, how 
 9069  many of the |εN |πrecords have passed before 
 9077  the sample is complete?)|'{A3}|5|5|≡6|≡.|9|4[|ε|*/M|↔p|↔M|\] 
 9082  |πWhat is the standard deviation of the value 
 9090  computed in the previous exercise?|'{A3}|5|5|≡7|≡.|9|4[|ε|*/M
 9095  |↔P|↔C|\] |πProve that any |εgiven |πchoice of 
 9102  |εn |πrecords from the set of |εN |πis obtained 
 9111  by Algorithm S with probability 1/(|ε|fN|d5n|)). 
 9117  |πTherefore the sample is completely unbiased.|'
 9123  |5|5|≡8|≡.|9|4[|ε|*/M|↔M|↔o|\] |πAlgorithm S computes 
 9127  one uniform deviate for each input record it 
 9135  handles. Find a more e∃cient way to determine 
 9143  the number of input records to skip before the 
 9152  _rst is selected, assuming that |εN/n |πis rather 
 9160  large. (We could iterate this process to select 
 9168  the remaining |εn|4α_↓|41 |πrecords.)|'|5|5|≡9|≡.|9|4[|ε|*/|↔
 9172  O|↔P|\] |πLet |εn|4α=↓|43. |πIf Algorithm R is 
 9179  applied to a _le containing 20 records numbered 
 9187  1 thru 20, and if the random numbers generated 
 9196  in step R3 are respectively|'{A6}4,|41,|46,|47,|45,|43,|45,|
 9201  411,|411,|43,|47,|49,|43,|411,|44,|45,|44,|;{A6}which 
 9203  records go into the reservoir? Which are in the 
 9212  _nal sample?|'{A3}|≡1|≡0|≡.|9|4[|ε|*/|↔O|↔C|\] 
 9215  |πModify Algorithm R so that the reservoir is 
 9223  eliminated, assuming that the |εn |πrecords of 
 9230  the current sample can be held in memory.|'{A3}|≡1|≡1|≡.|9|4
 9238  [|ε|*/M|↔P|↔C|\] |πLet |εp|βm |πbe the probability 
 9244  that exactly |εm |πelements are put into the 
 9252  reservoir during the _rst pass of Algorithm R. 
 9260  Determine the generating function |εG(z)|4α=↓|4|¬K|βmp|βmz|g
 9264  m, |πand _nd the mean and standard deviation. 
 9272  (Use the ideas of Section 1.2.10.)|'{A3}|≡1|≡2|≡.|9|4[|ε|*/M|
 9278  ↔P|↔o|\] |πThe gist of Algorithm P is that any 
 9287  permutation |ε|≤p |πcan be uniquely written as 
 9294  a product of transpositions in the form |ε|≤p|4α=↓|4(a|βtt)|
 9301  4.|4.|4.|4(a|β33)(a|β22), |πwhere |ε1|4|¬E|4a|βj|4|¬E|4j 
 9304  |πfor |εt|4|¬R|4j|4|¬Q|41. |πProve that there 
 9309  is also a unique representation of the form |ε|≤p|4α=↓|4(b|β
 9317  22)(b|β33)|4.|4.|4.|4(b|βtt), |πwhere |ε1|4|¬E|4b|βj|4|¬E|4j
 9319   |πfor |ε1|4|¬W|4j|4|¬E|4t, |πand give an algorithm 
 9326  which computes the |εb'|πs from the |εa'|πs in 
 9334  |εO(t) |πst{U0}{H10L12M29}W58320#Computer Programming|9(Knut
folio 190 galley 13
 9336  h/Addison-Wesley)|9f. 190|9ch. 3|9g. 13c|'{A6}{H9L11}|≡1|≡3|
 9340  ≡.|9|4[|ε|*/M|↔P|↔L|\] |π(S. W. Golomb.) A common 
 9346  procedure for card shu+ing is to divide the deck 
 9355  into two parts as equal as possible, and to ``ri+e'' 
 9365  them together. (See the discussion of card-playing 
 9372  etiquette in Hoyle's rules of card games; we 
 9380  read, ``A shu+e of this sort should be made about 
 9390  three times to mix the cards thoroughly.'') Consider 
 9398  a deck of |ε2n|4α_↓|41 |πcards |εX|β1,|4X|β2,|4.|4.|4.|4,|4X
 9403  |β2|βn|βα_↓|β1; |πa ``perfect shu+e'' |εs |πdivides 
 9409  this deck into |εX|β1,|4X|β2,|4.|4.|4.|4,|4X|βn 
 9413  |πand |εX|βn|βα+↓|β1,|4.|4.|4.|4,|4X|β2|βn|βα_↓|β1 
 9415  |πand perfectly interleaves them, to obtain |εX|β1, 
 9422  X|βn|βα+↓|β1, X|β2, X|βn|βα+↓|β2,|4.|4.|4.|4, 
 9425  X|β2|βn|βα_↓|β1,|4X|βn. |πThe ``cut'' operation 
 9429  |εc|gj |πchanges |εX|β1, X|β2,|4.|4.|4.|4, X|β2|βn|βα_↓|β1 
 9434  |πinto |εX|βj|βα+↓|β1,|4.|4.|4.|4, X|β2|βn|βα_↓|β1,|4X|β1,|4
 9436  .|4.|4.|4,|4X|βj. |πShow that by combining perfect 
 9442  shu+es and cuts, at most (|ε2n|4α_↓|41)(2n|4α_↓|42) 
 9448  |πdi=erent arrangements of the deck are possible, 
 9455  if |εn|4|¬Q|41.|'|≡1|≡4|≡.|9|4[|ε|*/|↔L|↔c|\] 
 9458  |π(Ole-Johan Dahl.) If |εX|βj|4α=↓|4j |πat the 
 9464  start of Algorithm P, and if we terminate that 
 9473  algorithm when |εj |πreaches the value |εt|4α_↓|4n, 
 9480  |πthe sequence |εX|βt|βα_↓|βn|βα+↓|β1,|4.|4.|4.|4,|4X|βt 
 9483  |πis a random permutation of a random combination 
 9491  of |εn |πelements. Show how to simulate the e=ect 
 9500  of this procedure using only |εO(n) |πcells of 
 9508  memory.|'{A24}{H10L12}|≤⊂|∨3|∨.|∨5|∨.|9|∨W|∨H|∨A|∨T|9|∨I|∨S|
 9509  9|∨A|9|∨R|∨A|∨N|∨D|∨O|∨M|9|∨S|∨E|∨Q|∨U|∨E|∨N|∨C|∨E|∨*3|'
 9510  {A12}|≡A|≡.|9|≡I|≡n|≡t|≡r|≡o|≡d|≡u|≡c|≡t|≡o|≡r|≡y 
 9511  |∨r|∨e|∨m|∨a|∨r|∨k|∨s|∨.|9|4We have seen in this 
 9516  chapter how to generate sequences|'{A9}|ε|¬DU|βn|¬F|4α=↓|4U|
 9521  β0,|4U|β1,|4U|β2,|4.|4.|4.|E|;(1)|?{A9}|πof real 
 9525  numbers in the range |ε0|4|¬E|4U|βn|4|¬W|41, 
 9530  |πand we have called them ``random'' sequences 
 9537  even though they are completely deterministic 
 9543  in character. To justify this terminology, we 
 9550  claimed that the numbers ``behave as if they 
 9558  are truly random.'' Such a statement may be satisfactory 
 9567  for practical purposes (at the present time), 
 9574  but it sidesteps a very important philosophical 
 9581  and theoretical question: Precisely what do we 
 9588  mean by ``random behavior''? A quantitative de_nition 
 9595  is needed. It is undesirable to talk about concepts 
 9604  that we do not really understand, especially 
 9611  since many apparently paradoxical statements 
 9616  can be made about random numbers.|'!|9|7The mathematical 
 9624  theory of probability and statistics carefully 
 9630  avoids answering the question; it refrains from 
 9637  making absolute statements, and instead expresses 
 9643  everything in terms of how much |εprobability 
 9650  |πis to be attached to statements involving random 
 9658  sequences of events. The axioms of probability 
 9665  theory are set up so that abstract probabilities 
 9673  can be computed readily, but nothing is said 
 9681  about what probability really signi_es, or how 
 9688  this concept can be applied meaningfully to the 
 9696  actual world. In the book |εProbability, Statistics, 
 9703  and Truth (|πMacmillan, 1957), R. von Mises discusses 
 9711  this situation in detail, and presents the view 
 9719  that a proper de_nition of probability depends 
 9726  on obtaining a proper de_nition of a rafdbmm 
 9734  s*?*?*?oper de_nition of a random sequence.|'!|9|7Let 
 9741  us paraphrase here some statements made by two 
 9749  of the many authors who have commented on the 
 9758  subject.|'{A12}{I1.7l}|εD. H. Lehmer (|*/|↔O|↔m|↔C|↔O|\)|*/: 
 9763  |\|π``A random sequence is a vague notion embodying 
 9771  the idea of a sequence in which each term is 
 9781  unpredictable to the uninitiated and whose digits 
 9788  pass a certain number of tests, traditional with 
 9796  statisticians and depending somewhat on the uses 
 9803  to which the sequence is to be put.``|'{A6}|εJ. 
 9812  N. Franklin (|*/|↔O|↔m|↔o|↔P|\)|*/: |\|π``The sequence 
 9817  (1) is random if it has every property that is 
 9827  shared by all in_nite sequences of independent 
 9834  samples of random variables from the uniform 
 9841  distribution.''|'{IC}{A12}!|9|7|πFranklin's statement 
 9844  essentially generalizes Lehmer's to say that 
 9850  the sequence must satisfy |εall |πstatistical 
 9856  tests. His de_nition is not completely precise, 
 9863  and we will see later that a reasonable interpretation 
 9872  of his statement leads us to conclude that there 
 9881  is no such thing as a random sequence*3 Let us 
 9891  begin with Lehmer's less restrictive statement 
 9897  and attempt to make |εit |πprecise. What we really 
 9906  want is a relatively short list of mathematical 
 9914  properties, each of which is satis_ed by our 
 9922  intuitive notion of a random sequence; furthermore, 
 9929  this list is to be complete enough so that we 
 9939  are willing to agree that |εany |πsequence satisfying 
 9947  these properties is ``random.'' In this section, 
 9954  we will develop what seems to be an adequate 
 9963  de_nition of randomness according to these criteria, 
 9970  although many interesting questions remain to 
 9976  be answered.|'!|9|7Let |εu |πand |εv |πbe real 
 9984  numbers, |ε0|4|¬E|4u|4|¬W|4v|4|¬E|41. |πIf |εU 
 9988  |πis a random variable which is uniformly distributed 
 9996  between 0 and 1, the probability that |εu|4|¬E|4U|4|¬W|4v 
10004  |πis equal to |εv|4α_↓|4u. |πFor example, the 
10011  probability that |f1|d33|)|4|¬E|4U|4|¬W|4|f2|d33|) 
10014  is |f1|d33|). |πHow can we translate this property 
10022  of the single number |εU |πinto a property of 
10031  the in_nite sequence |εU|β0,|4U|β1,|4U|β2,|4.|4.|4.|4? 
10035  |πThe obvious answer is to count how many times 
10044  |εU|βn |πlies between |εu |πand |εv, |πand the 
10052  average number of times should equal |εv|4α_↓|4u. 
10059  |πOur intuitive idea of probability is based 
10066  in this way on the frequency of occurrence. More 
10075  precisely, let |εv(n) |πbe the number of values 
10083  of |εj, 0|4|¬E|4j|4|¬W|4n, |πsuch that |εu|4|¬E|4U|βj|4|¬W|4
10088  v; |πwe want the ratio |ε|≤n(n)/n |πto approach 
10096  the value |εv|4α_↓|4u |πas |εn |πapproaches in_nity:|'
10103  {A9}|uwlim|uc|)|.|εn|¬M|¬X|)|4|≤n(n)/n|4α=↓|4v|4α_↓|4u.|;
10104  {A9}|πIf this condition holds for all choices 
10111  of |εu |πand |εv, |πthe sequence is said to be 
10121  |εequidistributed.|'!|9|7|πLet |εS(n) |πbe a 
10126  statement about the integer |εn |πand the sequence 
10134  |εU|β1,|4U|β2,|4.|4.|4.|4; |πfor example, |εS(n) 
10138  |πmight be the statement considered above, namely 
10145  ``u|4|¬E|4U|βn|4|¬W|4v.''  |πWe can generalize 
10150  the idea used in the preceding paragraph to de_ne 
10159  ``the probability that |εS(n) |πis true'' with 
10166  respect to a particular in_nite sequence: Let 
10173  |ε|≤n(n) |πbe the number of values of |εj, 0|4|¬E|4j|4|¬W|4n
10181  , |πsuch that |εS(j) |πis true.|'{A12}|≡D|≡e|≡_|≡n|≡i|≡t|≡i|
10187  ≡o|≡n |≡A|≡.|9|4|εWe say |πPr{H12}({H10}|εS(n){H12}){H10}|4α
10190  =↓|4|≤l, if |πlim|ε|βn|β|¬M|β|¬X|4|≤n(n)/n|4α=↓|4|≤l. 
10193  (|πRead, ``The probability that |εS(n) |πis true 
10200  is |ε|≤l, |πif the limit as |εn |πtends to in_nity 
10210  of |ε|≤n(n)/n |πis |ε|≤l.'')|'{A12}|πIn terms 
10216  of this notation, the sequence |εU|β0,|4U|β1,|4.|4.|4. 
10222  |πis equidistributed if and only if Pr(|εu|4|¬E|4U|βn|4|¬W|4
10228  v)|4α=↓|4v|4α_↓|4u, |πfor all real numbers |εu, 
10234  v |πwith |ε0|4|¬E|4u|4|¬W|4v|4|¬E|41.|'!|9|7|πA 
10238  sequence may be equidistributed without being 
10244  random. For example, if |εU|β0, U|β1,|4.|4.|4. 
10250  |πand |εV|β0,|4V|β1,|4.|4.|4. |πare equidistributed 
10254  sequences, it is not hard to show that the sequence|'
10264  {A9}|εW|β0,|4W|β1,|4W|β2,|4W|β3,|4.|4.|4.|4α=↓|4|f1|d32|)U|β
10264  0,|4|f1|d32|)|4α+↓|4|f1|d32|)V|β0,|4|f1|d32|)U|β1,|4|f1|d32|
10264  )|4α+↓|4|f1|d32|)V|β1,|4.|4.|4.|E|;(3)|?{A9}|πis 
10267  also equidistributed, since the sequence |ε|f1|d32|)U|βo,|4|
10272  f1|d32|)U|β1,|4.|4.|4. |πis equidistributed between 
10276  0 and |f1|d32|), |πwhile the alternate terms, 
10283  |f1|d32|)|4α+↓|4|f1|d32|)|εV|β0,|4|f1|d32|)|4α+↓|4|f1|d32|)V
10283  |β1,|4.|4.|4.|4, |πare equidistributed between 
10287  |f1|d32|) and 1. In the sequence of |εW|1'|πs, 
10295  a value less than |f1|d32|) is always followed 
10303  by a value greater than or equal to |f1|d32|), 
10312  and conversely; hence that sequence is not random 
10320  by any reasonable de_nition. A stronger property 
10327  than equidistribution is needed.|'!|9|7A natural 
10333  generalization of the equidistribution property, 
10338  which removes the objection stated in the preceding 
10346  paragraph, is to consider adjacent pairs of numbers 
10354  of our sequence. We can require the sequence 
10362  to satisfy the condition|'{A9}Pr(|εu|β1|4|¬E|4U|βn|4|¬W|4v|β
10366  1!!|πand!!|εu|β2|4|¬E|4U|βn|βα+↓|β1|4|¬W|4v|β2)|4α=↓|4(v|β1|
10366  4α_↓|4u|β1)(v|β2|4α_↓|4u|β2)|E|;(4)|?{A9}|πfor 
10369  any four numbers |εu|β1, v|β1, u|β2, v|β2 |πwith 
10377  |ε0|4|¬E|4u|β1|4|¬W|4v|β1|4|¬E|41,|40|4|¬E|4u|β2|4|¬W|4v|β2|
10377  4|¬E|41. |πIn general, for any positive integer 
10384  |εk |πwe can require our sequence to be |εk-distributed 
10393  |πin the following sense:|'{A12}|≡D|≡e|≡_|≡n|≡i|≡t|≡i|≡o|≡n 
10398  |≡B|≡.|9|4|εThe sequence (|*/|↔O|\) is said to 
10404  be k-distributed if|'{A9}|πPr(|εu|β1|4|¬E|4U|βn|4|¬W|4v|β1,|
10407  4.|4.|4.|4,|4u|βk|4|¬E{U0}{H10L12M29}W58320#Computer 
folio 194 galley 14
10408  Programming|9(Knuth/Addison-Wesley)|9f. 194|9ch. 
10410  3|9g. 14c|'{A6}!|9|7|πAn equidistributed sequence 
10415  is a 1-distributed sequence. Note that if |εk|4|¬Q|4, 
10423  |πa |εk-|πdistributed sequence is always (|εk|4α_↓|41)-|πdis
10428  tributed, since we may set |εu|βk|4α=↓|40 |πand 
10435  |εv|βk|4α=↓|41 |πin Eq. (5). Thus, in particular, 
10442  any sequence which is known to be 4-distributed 
10450  must also be 3-distributed, 2-distributed, and 
10456  equidistributed. We can investigate the largest 
10462  |εk |πfor which a given sequence is |εk-|πdistributed; 
10470  and this leads us to formulate|'{A12}|≡D|≡e|≡_|≡n|≡i|≡t|≡i|≡
10476  o|≡n |≡C|≡.|9|4|εA sequence is said to be |¬X-distributed 
10484  if it is k-distributed for all positive integers 
10492  k.|'{A12}!|9|7|πSo far we have considered ``[0,|41) 
10499  sequences,'' i.e., sequences of real numbers 
10505  lying between zero and one. The same ideas apply 
10514  to integer-valued sequences; let us say a sequence 
10522  |ε|¬DX|βn|¬F|4α=↓|4X|β0,|4X|β1,|4X|β2,|4.|4.|4. 
10523  |πis a ``|εb-|πary sequence'' if each |εX|βn 
10530  |πis one of the integers 0,|41,|4.|4.|4.|4,|4b|4α_↓|41. 
10536  |πThus, a 2-ary (binary) sequence is a sequence 
10544  of zeros and ones.|'!|9|7A |εk-|πdigit ``|εb-|πary 
10551  number'' |εx|β1x|β2|4.|4.|4.|4x|βk |πis an ordered 
10556  set of |εk |πintegers, where |ε0|4|¬E|4x|βj|4|¬W|4b 
10562  |πfor |ε1|4|¬E|4j|4|¬E|4k.|'{A12}|π|≡D|≡e|≡_|≡n|≡i|≡t|≡i|≡o|
10564  ≡n |≡D|≡.|9|4|εA b-ary sequence is said to be 
10572  k-distributed if|'{A9}|πPr(|εX|βnX|βn|βα+↓|β1|4.|4.|4.|4X|βn
10574  |βα+↓|βk|βα_↓|β1|4α=↓|4x|β1x|β2|4.|4.|4.|4x|βk)|4α=↓|41/b|gk
10574  |E|;(6)|?{A9}for all b-ary numbers x|β1x|β2|4.|4.|4.|4x|βk.|
10580  '{A12}!|9|7|πIt is clear from this de_nition 
10587  that if |εU|β0,|4U|β1,|4.|4.|4. |πis a |εk-|πdistributed 
10593  [0,|41) sequence, then the sequence |↔l|εbU|β0|↔L, 
10599  |↔lbU|β1|↔L,|4.|4.|4. |πis a |εk-|πdistributed 
10603  |εb-|πary sequence. (For if we set |εu|βj|4α=↓|4x|βj/b, 
10610  v|βj|4α=↓|4(x|βj|4α+↓|41)/b, X|βn|4α=↓|4|↔lbU|βn|↔L, 
10612  |πEq. (5) becomes Eq. (6).{H12}){H10} Furthermore, 
10618  every |εk-|πdistributed |εb-|πary sequence, for 
10623  |εk|4|¬Q|41, |πis also (|εk|4α_↓|41)-|πdistributed: 
10627  we add together the probabilities for the |εb-|πary 
10635  numbers |εx|β1|4.|4.|4.|4x|βk|βα_↓|β10, x|β1|4.|4.|4.|4x|βk|
10637  βα_↓|β11,|4.|4.|4.|4, x|β1|4.|4.|4.|4x|βk|βα_↓|β1(b|4α_↓|41)
10638   |πto obtain|'{A9}|πPr(|εX|βn|4.|4.|4.|4X|βn|βα+↓|βk|βα_↓|β2
10641  |4α=↓|4x|β1|4.|4.|4.|4x|βk|βα_↓|β1)|4α=↓|41/b|gk|gα_↓|g1.|;
10642  {A9}(|πProbabilities for disjoint events are 
10647  additive; see exercise 5.) It therefore is natural 
10655  to speak of an |¬X-distributed |εb-|πary sequence, 
10662  as in De_nition C above.|'!|9|7|πThe representation 
10669  of a positive real number in the radix-|εb |πnumber 
10678  system may be regarded as a |εb-|πary sequence; 
10686  for example, |ε|≤p |πcorresponds to the 10-ary 
10693  sequence 3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 
10706  9,|4.|4.|4.|4. |πIt has been conjectured that 
10712  this sequence is |¬X-distributed, but nobody 
10718  has yet been able to prove that it is even 1-distributed.|'
10729  !|9|7|πLet us analyze these concepts a little 
10736  more closely in the case when |εk |πequals a 
10745  million. A binary sequence which is 1000000-distributed 
10752  is going to have runs of a million zeros in a 
10763  row*3 Similarly, a [0,|41) sequence which is 1000000-distribu
10770  ted is going to have runs of a million consecutive 
10780  values each of which is less than |f1|d32|). 
10788  It is true that this will happen only (|f1|d32|))|g1|g0|g0|g
10796  0|g0|g0|g0 of the time, on the average, but the 
10805  fact is that it |εdoes |πhappen. Indeed, this 
10813  phenomenon will occur in any truly random sequence, 
10821  using our intuitive notion of ``truly random.'' 
10828  One can easily imagine that such a situation 
10836  will have a drastic e=ect if this set of a million 
10847  ``truly random'' numbers is being used in a computer-simulat
10855  ion experiment; there would be good reason to 
10863  complain about the random-number generator*3 However, 
10869  if we have a sequence of numbers which never 
10878  has runs of a million consecutive |εU'|πs less 
10886  than |f1|d32|), the sequence is not random, and 
10894  it will not be a suitable source of numbers for 
10904  other conceivable applications which use extremely 
10910  long blocks of |εU'|πs as input. In summary, 
10918  |εa truly random sequence will exhibit local 
10925  nonrandomness|*/; |\|πlocal nonrandomness is necessary 
10930  in some applications, but it is disastrous in 
10938  others. We are forced to conclude that |εno sequence 
10947  of ``random'' numbers can be adequate for every 
10955  application.|'!|9|7|πIn a similar vein, one may 
10962  argue that there is no way to judge whether a 
10972  |ε⊂nite |πsequence is random or not; any particular 
10980  sequence is just as likely as any other one. 
10989  These facts are de_nitely stumbling blocks if 
10996  we are ever to have a useful de_nition of randomness, 
11006  but they are not really cause for alarm. It is 
11016  still possible to give a de_nition for the randomness 
11025  of in_nite sequences of real numbers in such 
11033  a way that the corresponding theory (viewed properly) 
11041  will give us a great deal of insight concerning 
11050  the ordinary _nite sequences of rational numbers 
11057  which are actually generated on a computer. Furthermore, 
11065  we shall see later in this section that there 
11074  are several plausible de_nitions of randomness 
11080  for _nite sequences are in fact available.|'{A12}|≡B|≡.|9|¬X
11087  |≡-|≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡e|≡d |≡s|≡e|≡q|≡u|≡e|≡n|≡c|≡
11088  e|≡s|≡.|9|4Let us now undertake a brief study 
11095  of the theory of |¬X-distributed sequences. To 
11102  describe the theory adequately, we will need 
11109  to use a little ``higher mathematics,'' so we 
11117  assume in the remainder of this subsection that 
11125  the reader knows the material ordinarily taught 
11132  in an ``advanced calculus'' course.|'|Ha*?*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!
11137